Vol de fusée : Schéma d'interpolation et discrétisation

Nous continuons le problème du vol de la fusée en abordant les schémas d'interpolation. Le choix du schéma d'interpolation est, à mon avis, une des deux étapes les plus importantes d'une résolution numérique, l'autre étant la discrétisation en elle-même. Le schéma choisi a une influence directe sur la précision, la vitesse de convergence, la stabilité et sur la difficulté d'implémentation.

Avant d'entrer dans le vif du sujet, on va établir quel est l'objectif d'un schéma d'interpolation. De façon générale, un schéma d'interpolation est une équation algébrique qu'on utilise pour représenter la variation d'une variable entre deux valeurs discrètes. Dans le cas des méthodes numériques, les schémas d'interpolation sont utilisés à deux reprises. Premièrement, ils servent à évaluer les dérivées dans les systèmes d'équations partielles afin de transformer ces systèmes en équations algébriques. Deuxièmement, ils peuvent aussi être utilisés afin de produire des visualisation de la solution. Pour l'instant, c'est la première utilité qui nous intéresse.

Dans le cas du vol de fusée, les seules dérivées dans les équations différentielles sont des dérivées temporelles. C'est donc dire que la seule interpolation que nous aurons à faire sera temporelle. Nous considérerons d'abord le cas où on discrétise le temps en intervalles constants. Nous appelerons cet intervalle $\Delta t$. Nous utiliserons une notation indicielle pour les variables discrètes.

$$t_n=n \Delta t$$ $$v_n=v(t_n)$$ $$h_n=h(t_n)$$

La manière la plus simple de décrire la variation d'une variable entre deux points, est de considérer une variation linéarire. Si on considère un nombre $a$ entre 0 et 1, la valeur de $t$ entre $t_n$ et $t_{n+1}$ peut s'écrire de la façon suivante.

$$t= \left( 1-a \right) t_n + a t_{n+1} = t_n + a \Delta t $$

Si on dérive cette expression par rapport à $t$, on obtient:

$$\frac{\mathrm{d} a}{\mathrm{d} t}=\frac{1}{\Delta t}$$

On considère la même variation linéaire entre $v_n$ et $v_{n+1}$.

$$v= \left( 1-a \right) v_n + a v_{n+1}$$

On veut maintenant obtenir la dérivée de $v$.

$$\frac{\mathrm{d} v}{\mathrm{d} t}=-\frac{\mathrm{d} a}{\mathrm{d} t} v_n+\frac{\mathrm{d} a}{\mathrm{d} t} v_{n+1}$$ $$\frac{\mathrm{d} v}{\mathrm{d} t}=\frac{v_{n+1}-v_n}{\Delta t}$$

On a développer une expression pour la dérivée d'une variable entre deux points discrèts. On peut substituter cette expression dans l'équation différentielle de vitesse.

$$\frac{v_{n+1}-v_n}{\Delta t}=-9.81+\frac{6500}{150-20t}-\frac{3 \cdot 1.091 \pi}{160 \left( 150 - 20t \right)} v^2$$ On remarque maintenant qu'il faut choisir une valeur discrète de $t$ et de $v$ pour le terme de résistance de l'air et le terme de propulsion. On peut utiliser l'interpolation linéaire, en fonction des valeurs discrètes, et du paramètre $a$. $$\frac{v_{n+1}-v_n}{\Delta t}=-9.81+\frac{6500}{ 150 - 20 \left( t_n + a \Delta t \right) }-\frac{3 \cdot 1.091 \pi}{160 \Big( 150 - 20 \left( t_n + a \Delta t \right) \Big)} {\Big( \left( 1-a \right) v_n + a v_{n+1} \Big)}^2$$

Le choix du paramètre correspond à choisir l'instant utilisé pour évaluer la résistance de l'air et la poussée. La valeur du paramètre $a$ qu'on choisi détermine si on utilise un schéma tmeporel explicite ($a=0$), implicite ($a=1$) ou semi-implicite ($ 0 < a < 1 $). Le cas $a=\frac{1}{2}$ correspond au schéma de Crank-Nicolson.

Ce qu'on veut faire par la suite avec l'expression est d'isoler $t_{n+1}$ afin d'avoir une solution qui progresse dans le temps (time-marching solution). Puisqu'on connait la valeur des variables à $t_0$, on pourra d'abord calculer avec notre formule la valeur des variable à $t_1$, et successivement. Il est donc algébriquement plus simple d'utiliser $a = 0$, puisque ce choix élimine le terme en $v_{n+1}$ dans le membre de droite.

$$\frac{v_{n+1}-v_n}{\Delta t}=-9.81+\frac{6500}{150-20t_n}-\frac{3 \cdot 1.091 \pi}{160 \left( 150 - 20t_n \right)} v_n^2$$ $$v_{n+1}=v_n + \Delta t \Big( -9.81+\frac{6500}{150-20t_n}-\frac{3 \cdot 1.091 \pi}{160 \left( 150 - 20t_n \right)} v_n^2 \Big) $$

Nous avons développé une expression algébrique qui permet de calculer la valeur de $v$ par incrément de $\Delta t$. Un schéma temporel linéaire explicite à été utilisé. On peut effectuer le même travail pour les autres équations différentielles de $v$ et pour l'équation différentielles de $h$.

$$v_{n+1}=v_n + \Delta t \Big( -9.81 -\frac{3 \cdot 1.091 \pi}{8000} v_n^2 \Big) $$ $$v_{n+1}=v_n + \Delta t \Big( -9.81 +\frac{3 \cdot 1.091 \pi}{8000} v_n^2 \Big) $$ $$h_{n+1}=h_n + \Delta t v_n $$

Les formules développées permettront de programmer facilement un algorithme de résolution numérique. En fonction de la performance de cet algorithme, d'autres schémas d'interpolation pourraient être implémentés.

Vol de fusée : Description des équations

Dans le billet précédent, nous avons commencé le problème du vol de fusée, tel que présenté dans le cours MAE6286 Practical Numerical Methods with Python. La définition du problème se trouve ici. Les fonctions pour la masse et pour le débit massique de carburant ont été déterminées. Afin de représenter le comportement de ces variables au moment où la fusée est à cours de carburant, ces fonctions sont définies par morceaux. Les vitesses terminales en ascension et en descente ont été également calculées. Ces vitesses correspondent aux limites que la fusée ne peut pas dépasser. La fusée n'atteindra pas nécessairement ces vitesses dans notre modèle numérique, cependant, un problème sera évident si la vitesse de la fusée est en dehors de l'intervalle que nous avons défini.

Dans ce billet, on écrira les équations différentielles à résoudre pour la vitesse et l'altitude. On note d'abord que les équations différentielles ne sont pas couplées entres elles. On pourra d'abord résoudre l'équation pour la vitesse, et ensuite, résoudre celle de l'altitude.

Puisque certaines variables de l'équation différentielle de la vitesse sont définies par morceaux, on devra résoudre l'équation de la vitesse par morceaux. C'est-à-dire qu'on aura une équation différentielles pour chaque morceau.

Premièrement, pour les 5 premières secondes, on peut substituer les variables de masse et de débit massique de carburant. On sait aussi que la vitesse est positive, ce qui nous permet de traiter la valeur absolue. On substitue aussi les valeurs numériques. $$\left (150-20t \right )\frac{\mathrm{d} v}{\mathrm{d} x}=-\left ( 150-20t \right )\cdot 9.81+20\cdot 325-\frac{1}{2}\cdot 1.091 \cdot \pi \cdot 0.5^2 \cdot 0.15 \cdot v^2$$ Finalement, on isole la dérivée de la vitesse. $$\frac{\mathrm{d} v}{\mathrm{d} t}=-9.81+\frac{6500}{150-20t}-\frac{3 \cdot 1.091 \pi}{160 \left( 150 - 20t \right)} v^2$$ Le premier terme du membre de droite correspond a l'accélération gravitationnelle. Le second correspond à l'accélération due à la poussée de la propulsion. Le dernier terme, non-linéaire, correspond à la résistance de l'air. Une solution exacte serait relativement facile à obtenir sans la présence de ce troisième terme.

Lorsque la solution de cette équation sera résolue pour les cinq premières secondes, l'état des variables correspondra à une seconde équation différentielle. Afin de traiter la valeur absolue, on considéra cette équation pendant que la vitesse est positive. Au premier instant où la vitesse est négative, on devra changer d'équation. $$\frac{\mathrm{d} v}{\mathrm{d} x}=-9.81-\frac{3 \cdot 1.091 \pi}{8000} v^2$$ L'équation différentielle a subi quelques changements, mais est encore relativement semblable. La masse ne dépend plus du temps, et le terme de poussée a disparu. On note ici que l'accélération doit être négative, et donc qu'un régime permanent n'est pas possible.

Lorsque la vitesse aura atteint son altitude maximale : c'est-à-dire lorsque sa vitesse commencera à être négative, il sera nécessaire de résoudre une troisième équation. $$\frac{\mathrm{d} v}{\mathrm{d} x}=-9.81+\frac{3 \cdot 1.091 \pi}{8000} v^2$$ On remarque ici que seul le signe de la résistance de l'air a changé.

En parallèle de l'équation de vitesse, on doit aussi résoudre l'équation de l'altitude. $$\frac{\mathrm{d} h}{\mathrm{d} x}=v$$ Le problème se termine au moment de l'écrasement, c'est-à-dire quand l'altitude est nulle est que la vitesse est négative. Dans ce billet, nous avons décrit les différentes équations différentielles que nous allons résoudre numériquement. Dans le prochain billet, nous aborderons la discrétisation de ces équations.

Vol de fusée : un premier regard

Le premier problème soumis durant le cours MAE6286 Practical Numerical Methods with Python est le problème du vol de fusée. Il est possible de consulter la définition du problème tel que présenté dans le cours en suivant ce lien.

Premièrement, la consommation de carburant est découplée du système d'équations différentielles. On peut donc traiter cet aspect du problème de façon indépendante. On peut écrire une fonction définie par morceaux pour le débit de carburant, en fonction du temps. $$\dot{m}_{p}(t) = \bigg \{ \begin{matrix} 20 & si & 0\leq t < 5\\ 0 & si & t \geq 5 \end{matrix}$$ À l'aide de cette fonction, on peut calculer la masse de carburant dans la fusée en fonction du temps, par intégration. $$m_{p}(t) = \bigg \{ \begin{matrix} 100-20t & si & 0\leq t < 5\\ 0 & si & t \geq 5 \end{matrix}$$ On peut donc avoir, à tout instant la masse d'essence en fonction du temps. On peut, entre autres, calculer la masse de carburant à 3.6 secondes.

Deuxièmement, on peut calculer la vitesse termimale de la fusée lors de son ascension et lors de sa descente. La connaissance de ces vitesses limites pourra servir à vérifier le fonctionnement de notre modèle numérique : en effet si ces vitesses sont dépassées, on saura qu'il y a un problème. Évidemment, le fait de ne pas dépasser ces vitesses ne veut cependant pas dire que notre modèle est bon...

Nous abordons d'abord l'ascencion. Par définition de la vitesse terminale, lors de l'ascension, on a deux relations : $$\begin{matrix} \frac{\mathrm{d} v}{\mathrm{d} x}=0 & et & \left | v \right |>0 \end{matrix}$$ L'équation différentielle de la vitesse devient une simple expression algébrique. $$\frac{1}{2}\rho C_D A v_{max}^{2} =\dot{m}_p v_e-(m_s+m_p)g$$ Selon cette équation, on constate que la vitesse maximale est la plus grande si le débit de carburant est maximal et si la masse de carburant est minimale. Or, c'est deux situations existent simultanément à t = 5 secondes. On peut remplacer par les valeurs numériques est résoudre pour v. $$v_{max} = 305.8 m/s$$ De facon similaire, lorsque lors de la descente, on a : $$ \begin{matrix} \frac{\mathrm{d} v}{\mathrm{d} x}=0 & , & \left | v \right |<0 & , & \dot{m}_p=0 & et & m_p=0 \end{matrix} $$ On obtient alors : $$\frac{1}{2}\rho C_D A v_{min}^{2} =m_s g$$ $$v_{min} = -87.4 m/s$$ Dans ce premier billet sur le problème du vol de fusée, on a décrit le débit massique et la masse du carburant en fonction du temps, à l'aide d'une fonction par morceau. Nous avons aussi établi les bornes inférieures et supérieures de la vitesse de la fusée.

Introduction

J'ouvre aujourd'hui un blogue. C'est essentiellement ma première expérience de publication, en général, et sur internet, en particulier. Je vais profiter de ce premier billet pour me présenter brièvement et donner un aperçu du contenu que je compte publier sur ce blogue.

Je m'appelle Alexandre Caron et habite à Montréal. J'ai complété en 2013 un Baccalauréat en génie mécanique à l'École de technologie supérieure et je suis présentement candidat à la maitrise en génie mécanique au même établissement d'enseignement. À titre d'information, mon projet de maitrise porte sur l'automatisation des calculs de stabilité aéroélastique, pour des raisons de confidentialité, le sujet de ma maitrise ne sera pas abordé directement dans ce blogue.

Au courant de mes études, j'ai toujours eu un intérêt pour les mathématiques. Pour la petite histoire, et pour rendre hommage aux individus concernés, il m'apparait évident que cet intérêt est largement attribuable aux différents enseignants de mathématiques que j'ai eu à travers mon parcours scolaire. Je me permet de nommer messieurs Albert Beaulieu et Olivier Arsenault qui ont particulièrement contribué au développement de ma curiosité pour les mathématiques. Le second est également responsable de ma découverte de la programmation, ce qui a certainement été un point tournant de mon approche des problèmes de natures mathématiques.

Dans ce blogue, je prévois publier du contenu portant sur l'analyse numérique des systèmes physiques. Pour ce faire, je compte présenter des exemples que je réalise, mais aussi des courts résumés de théories.

Mon objectif, en tenant un blogue, est de créer un ensemble de notes ordonnées et claires portant sur les différents sujets que j'aborde. En réalité, je crée ce contenu pour moi-même, et tant mieux si cela peut être utile et intéressant pour d'autres.

L'élément déclencheur qui m'a emmené à démarrer un blogue est mon inscription au cours en ligne MAE6286 Practical Numerical Methods with Python. C'est un cours qui utilise le format ''Formation en ligne ouverte à tous'' et qui est offert par un groupe de professeurs qui a à cœur de transmettre leurs connaissances au plus grand nombre. J'ai préalablement découvert de cours à travers le blogue du professeur Lorena A, Barba et sa série 12 Steps to Navier-Stokes. Je recommande d'ailleurs fortement cette série, qui présente très bien un exemple d'application du calcul numérique.

À travers les différentes publications, j'emploierai certainement plusieurs langages de programmation. Le cours MAE6286 est présenté à l'aide de Python, j'utiliserai cette opportunité pour apprendre les rudiments de ce langage, dans son application pour le calcul numérique. J'ai déjà réalisé des projets assez complexes avec Matlab et Fortran. Finalement, je possède des notions de C et de C++ assez limitées, mais que j'ai l'intention de développer. J'espère présenter des billets en utilisant ces différents langages.

Je vous souhaite à tous la bienvenue sur mon blogue, j'espère que vous vous y plairez.

Au plaisir de discuter avec vous,

Alexandre