Schéma monotone décentré amont

Summarize Timeline Top Qs Fact Check

On étudie les principes fondamentaux du schéma MUSCL en considérant le système scalaire 1D simple du premier ordre suivant, dans lequel on suppose qu'une onde se propage dans la direction positive :

${\displaystyle u_{t}+F_{x}\left(u\right)=0}$

Où ${\displaystyle u}$ représente une variable d'état et ${\displaystyle \textstyle F}$ représente son flux.

Le schéma de Godounov utilise des approximations constantes par morceaux pour chaque cellule et aboutit à une discrétisation décentrée du premier ordre du problème ci-dessus, les centres des cellules étant indexés par ${\displaystyle i}$. Un schéma semi-discret peut être défini comme suit :

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}+{\frac {1}{\Delta x_{i}}}\left[F\left(u_{i}\right)-F\left(u_{i-1}\right)\right]=0}$

Ce schéma de base ne permet pas de traiter les chocs ni les discontinuités abruptes, car ces dernières ont tendance à s'estomper. Un exemple de cet effet est illustré dans la figure ci-contre, qui représente une équation d'advection unidimensionnelle avec une onde en échelon se propageant vers la droite. La simulation a été réalisée avec un maillage de 200 cellules et un intégrateur temporel de Runge-Kutta d'ordre 4 (RK4).

Pour une meilleure résolution des discontinuités, le schéma de Godunov peut être étendu en utilisant des approximations linéaires par morceaux de chaque cellule, ce qui conduit à un schéma aux différences finies centrées d'ordre 2 en espace. Les approximations linéaires par morceaux sont obtenues à partir de :

${\displaystyle u\left(x\right)=u_{i}+{\frac {\left(x-x_{i}\right)}{\left(x_{i+1}-x_{i}\right)}}\left(u_{i+1}-u_{i}\right)\qquad \forall x\in [x_{i},x_{i+1}]}$

Ainsi, en évaluant les flux aux bords de la cellule, on obtient le schéma semi-discret suivant :

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}+{\frac {1}{\Delta x_{i}}}\left[F\left(u_{i+1/2}\right)-F\left(u_{i-1/2}\right)\right]=0}$

où ${\displaystyle u_{i+1/2}}$ et ${\displaystyle u_{i-1/2}}$ sont les valeurs approchées par morceaux des variables de bord de cellule :

${\displaystyle u_{i+1/2}=1/2\left(u_{i}+u_{i+1}\right)}$

${\displaystyle u_{i-1/2}=1/2\left(u_{i-1}+u_{i}\right)}$

Bien que le schéma du second ordre présenté ci-dessus offre une meilleure précision pour les solutions régulières, il ne s'agit pas d'un schéma TVD et il introduit des oscillations parasites dans la solution en présence de discontinuités ou de chocs. Un exemple de cet effet est illustré dans la figure ci-contre, qui représente une équation d'advection unidimensionnelle ${\displaystyle \textstyle u_{t}+u_{x}=0}$ avec une onde en échelon se propageant vers la droite. Cette perte de précision est prévisible en vertu du théorème de Godounov. La simulation a été réalisée avec un maillage de 200 cellules et a utilisé RK4 pour l'intégration temporelle.

Les schémas numériques basés sur MUSCL étendent l'idée d'utiliser une approximation linéaire par morceaux pour chaque cellule en utilisant des états extrapolés à gauche et à droite en utilisant un limiteur de pente. Il en résulte le schéma de discrétisation TVD haute résolution suivant :

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}+{\frac {1}{\Delta x_{i}}}\left[F\left(u_{i+1/2}^{*}\right)-F\left(u_{i-1/2}^{*}\right)\right]=0}$

Ce schéma peut également s'écrire sous la forme plus concise suivante :

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}+{\frac {1}{\Delta x_{i}}}\left[F_{i+1/2}^{*}-F_{i-1/2}^{*}\right]=0}$

Les flux numériques ${\displaystyle F_{i\pm 1/2}^{*}}$ correspondent à une combinaison non linéaire d'approximations du premier et du second ordre de la fonction de flux continue.

Les symboles ${\displaystyle u_{i+1/2}^{*}}$ et ${\displaystyle u_{i-1/2}^{*}}$ représentent des fonctions dépendant du schéma (des variables extrapolées limitées aux bords de la cellule), c'est-à-dire (indice L pour left, R pour right):

${\displaystyle u_{i+1/2}^{*}=u_{i+1/2}^{*}\left(u_{i+1/2}^{L},u_{i+1/2}^{R}\right)\,,\;u_{i-1/2}^{*}=u_{i-1/2}^{*}\left(u_{i-1/2}^{L},u_{i-1/2}^{R}\right),}$

où, en utilisant les pentes sous le vent :

${\displaystyle u_{i+1/2}^{L}=u_{i}+1/2\,\phi \left(r_{i}\right)\left(u_{i+1}-u_{i}\right)\,,\;u_{i+1/2}^{R}=u_{i+1}-1/2\,\phi \left(r_{i+1}\right)\left(u_{i+2}-u_{i+1}\right)}$

${\displaystyle u_{i-1/2}^{L}=u_{i-1}+1/2\,\phi \left(r_{i-1}\right)\left(u_{i}-u_{i-1}\right)\,,\;u_{i-1/2}^{R}=u_{i}-1/2\,\phi \left(r_{i}\right)\left(u_{i+1}-u_{i}\right),}$

et

${\displaystyle r_{i}=(u_{i}-u_{i-1})/(u_{i+1}-u_{i})}$

La fonction ${\displaystyle \textstyle \phi (r_{i})}$ est un limiteur de flux qui assure la discrétisation TVD de la solution, évitant ainsi les oscillations parasites qui sans cela se produiraient autour des discontinuités ou des chocs. Le limiteur est nul pour ${\displaystyle \textstyle r\leq 0}$ et égal à 1 pour ${\displaystyle \textstyle r=1}$. Ainsi, la précision d'une discrétisation TVD se dégrade au premier ordre aux extrema locaux, mais tend vers le second ordre sur les parties régulières du domaine.

L'algorithme est simple à implémenter. Une fois un schéma approprié pour ${\displaystyle \textstyle F_{i+1/2}^{*}}$ choisi, tel que le schéma de Kurganov et Tadmor (voir ci-dessous), la résolution peut être effectuée à l'aide de techniques d'intégration numérique classiques.

Summarize Timeline Top Qs Fact Check

Un précurseur du schéma centré de Kurganov et Tadmor (KT)^[4] est le schéma centré décalé de Nessyahu et Tadmor (NT)^[5]. Il s'agit d'un schéma du second ordre, sans solveur de Riemann, à haute résolution, utilisant la reconstruction MUSCL. C'est une méthode entièrement discrète, simple à implémenter, applicable aux problèmes scalaires et vectoriels. Elle peut être vue comme un flux de Rusanov (également appelé flux de Lax-Friedrichs local) complété par des reconstructions d'ordre élevé. L'algorithme est basé sur la méthode des différences finies ou des différences centrées et offre des performances comparables aux solveurs de type Riemann pour la résolution d'équations aux dérivées partielles décrivant des systèmes présentant des phénomènes à fort gradient.

Le schéma KT étend le schéma NT et présente une viscosité numérique plus faible que ce dernier. Il offre également l'avantage de pouvoir être implémenté en mode discret ou semi-discret. Nous considérons ici le schéma semi-discret.

Le calcul des flux est le suivant :

${\displaystyle F_{i-1/2}^{*}=1/2\left\{\left[F\left(u_{i-1/2}^{R}\right)+F\left(u_{i-1/2}^{L}\right)\right]-a_{i-1/2}\left[u_{i-1/2}^{R}-u_{i-1/2}^{L}\right]\right\}}$

${\displaystyle F_{i+1/2}^{*}=1/2\left\{\left[F\left(u_{i+1/2}^{R}\right)+F\left(u_{i+1/2}^{L}\right)\right]-a_{i+1/2}\left[u_{i+1/2}^{R}-u_{i+1/2}^{L}\right]\right\}}$

où la vitesse de propagation locale ${\displaystyle \textstyle a_{i\pm 1/2}}$ est la valeur absolue maximale de la valeur propre du jacobien de ${\displaystyle \textstyle F\left(u\left(x,t\right)\right)}$ sur les cellules ${\displaystyle ({i},{i\pm 1})}$ donnée par :

${\displaystyle a_{i+1/2}\left(t\right)=\max \left[\rho \left({\frac {\partial F\left(u_{i+1/2}^{L}\left(t\right)\right)}{\partial u}}\right)\,,\;\rho \left({\frac {\partial F\left(u_{i+1/2}^{R}\left(t\right)\right)}{\partial u}}\right)\right]}$

et ${\displaystyle \rho \left({\frac {\partial F\left(u\left(t\right)\right)}{\partial u}}\right)}$ représente le rayon spectral de ${\displaystyle {\frac {\partial F\left(u\left(t\right)\right)}{\partial u}}}$.

Au-delà des vitesses liées à la condition de Courant-Friedrichs-Lewy (CFL), aucune information caractéristique n'est requise.

Le calcul de flux présenté ci-dessus est généralement appelé « méthode de Lax-Friedrichs » Il convient toutefois de noter que cette expression de flux n'apparaît pas dans Lax (1954), mais plutôt dans Rusanov (1961).

L'efficacité d'un schéma à haute résolution est illustrée par le diagramme ci-contre, qui représente l'équation d'advection unidimensionnelle ${\displaystyle \textstyle u_{t}+u_{x}=0}$ avec une onde en échelon se propageant vers la droite. La simulation a été réalisée sur un maillage de 200 cellules, en utilisant le schéma central de Kurganov et Tadmor avec limiteur de flux Superbee et en employant RK4 pour l'intégration temporelle. Ce résultat de simulation contraste fortement avec les résultats obtenus précédemment par la méthode des différences finies décentrées d'ordre 1 et par la méthode des différences finies centrées d'ordre 2. Ce schéma donne également de bons résultats lorsqu'il est appliqué à des systèmes d'équations ; voir ci-dessous les résultats obtenus avec ce schéma appliqué aux équations d'Euler. Toutefois, il convient de choisir avec soin un limiteur approprié car, par exemple, le limiteur Superbee peut engendrer une variation excessivement rapide pour certaines ondes lisses.

Le schéma peut facilement intégrer des termes de diffusion, s'ils sont présents. Par exemple, si l'on étend le problème scalaire 1D ci-dessus pour inclure un terme de diffusion, on obtient :

${\displaystyle u_{t}+F_{x}\left(u\right)=Q_{x}\left(u,u_{x}\right)}$

pour lequel Kurganov et Tadmor proposent l'approximation par différences centrées suivante :

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}=-{\frac {1}{\Delta x_{i}}}\left[F_{i+1/2}^{*}-F_{i-1/2}^{*}\right]+{\frac {1}{\Delta x_{i}}}\left[P_{i+1/2}-P_{i-1/2}\right]}$

où

${\displaystyle P_{i+1/2}=1/2\left[Q\left(u_{i},{\frac {u_{i+1}-u_{i}}{\Delta x_{i}}}\right)+Q\left(u_{i+1},{\frac {u_{i+1}-u_{i}}{\Delta x_{i}}}\right)\right]}$

${\displaystyle P_{i-1/2}=1/2\left[Q\left(u_{i-1},{\frac {u_{i}-u_{i-1}}{\Delta x_{i-1}}}\right)+Q\left(u_{i},{\frac {u_{i}-u_{i-1}}{\Delta x_{i-1}}}\right)\right]}$

Les détails complets de l'algorithme (versions « complète » et « semi-discrète ») et sa dérivation sont disponibles dans l'article original, ainsi que plusieurs exemples 1D et 2D^[4]. Des informations complémentaires sont également disponibles dans l'article antérieur de Nessyahu et Tadmor^[5].

Remarque : Ce schéma a été initialement présenté par Kurganov et Tadmor comme un schéma du second ordre basé sur une « extrapolation linéaire ». Un article ultérieur^[6] démontre qu'il peut également servir de base à un schéma du troisième ordre. Un exemple d'advection 1D et un exemple d'équation d'Euler de leur schéma, utilisant une reconstruction parabolique (3e ordre), sont présentés dans les sections « reconstruction parabolique » et « équation d'Euler » ci-dessous.

Summarize Timeline Top Qs Fact Check

Il est possible d'étendre l'idée d'extrapolation linéaire à une reconstruction d'ordre supérieur, comme illustré dans le diagramme ci-contre. Dans ce cas, les états gauche et droit sont estimés par interpolation d'une équation aux différences finies décentrée d'ordre deux. On obtient ainsi un schéma de reconstruction parabolique d'ordre trois en espace.

En suivant l'approche de Kermani^[7] on présente un schéma décentrée d'ordre trois, où les symboles ${\displaystyle \textstyle u_{i+1/2}^{*}}$ et ${\displaystyle \textstyle u_{i-1/2}^{*}}$ représentent des fonctions dépendant du schéma (des variables de bord de cellule reconstruites, de taille limitée). Mais dans ce cas, elles reposent sur des états reconstruits paraboliquement, c'est-à-dire :

${\displaystyle u_{i+1/2}^{*}=f\left(u_{i+1/2}^{L},u_{i+1/2}^{R}\right),\quad u_{i-1/2}^{*}=f\left(u_{i-1/2}^{L},u_{i-1/2}^{R}\right)}$

et

${\displaystyle u_{i+1/2}^{L}=u_{i}+{\frac {\phi \left(r_{i}\right)}{4}}\left[\left(1-\kappa \right)\delta u_{i-1/2}+\left(1+\kappa \right)\delta u_{i+1/2}\right]}$

${\displaystyle u_{i+1/2}^{R}=u_{i+1}-{\frac {\phi \left(r_{i+1}\right)}{4}}\left[\left(1-\kappa \right)\delta u_{i+3/2}+\left(1+\kappa \right)\delta u_{i+1/2}\right]}$

${\displaystyle u_{i-1/2}^{L}=u_{i-1}+{\frac {\phi \left(r_{i-1}\right)}{4}}\left[\left(1-\kappa \right)\delta u_{i-3/2}+\left(1+\kappa \right)\delta u_{i-1/2}\right]}$

${\displaystyle u_{i-1/2}^{R}=u_{i}-{\frac {\phi \left(r_{i}\right)}{4}}\left[\left(1-\kappa \right)\delta u_{i+1/2}+\left(1+\kappa \right)\delta u_{i-1/2}\right]}$

où ${\displaystyle \textstyle \kappa =1/3}$ et :

${\displaystyle \delta u_{i+1/2}=u_{i+1}-u_{i}}$

${\displaystyle \delta u_{i-1/2}=u_{i}-u_{i-1}}$

${\displaystyle \delta u_{i+3/2}=u_{i+2}-u_{i+1}}$

${\displaystyle \delta u_{i-3/2}=u_{i-1}-u_{i-2}}$

et le limiteur ${\displaystyle \textstyle \phi \left(r\right)}$ est identique à celui présentée précédemment.

La reconstruction parabolique est simple à implémenter et peut être utilisée avec le schéma de Kurganov et Tadmor à la place de l'extrapolation linéaire mentionnée plus haut. Ceci a pour effet d'élever la solution spatiale du schéma KT au 3e ordre. Elle donne de bons résultats pour la résolution des équations d'Euler (voir ci-dessous). Cette augmentation de l'ordre spatial présente certains avantages par rapport aux schémas du 2e ordre pour les solutions régulières ; cependant, pour les chocs, elle est plus dissipative (comparer le diagramme ci-contre avec la solution ci-dessus obtenue à l'aide de l'algorithme KT avec extrapolation linéaire et limiteur Superbee). Cette simulation a été réalisée sur un maillage de 200 cellules en utilisant le même algorithme KT, mais avec une reconstruction parabolique. L'intégration temporelle a été effectuée par RK4 et la forme alternative du limiteur de van Albada ${\displaystyle \textstyle \phi _{va}(r)={\frac {2r}{1+r^{2}}}}$ a été utilisée pour éviter les oscillations parasites.

Summarize Timeline Top Qs Fact Check

Par souci de simplicité, nous considérons le cas unidimensionnel sans transfert de chaleur ni force volumique. Par conséquent, sous forme de vecteurs de conservation, les équations d'Euler se réduisent à :

${\displaystyle {\frac {\partial \mathbf {U} }{\partial t}}+{\frac {\partial \mathbf {F} }{\partial x}}=0}$

où

${\displaystyle \mathbf {U} ={\begin{pmatrix}\rho \\\rho u\\E\end{pmatrix}}\qquad \mathbf {F} ={\begin{pmatrix}\rho u\\p+\rho u^{2}\\u(E+p)\end{pmatrix}},\qquad }$

et où ${\displaystyle \textstyle \mathbf {U} }$ est un vecteur d'état et ${\displaystyle \textstyle \mathbf {F} }$ est un vecteur de flux.

Les équations ci-dessus représentent la conservation de la masse, de la quantité de mouvement et de l'énergie. Il y a donc trois équations et quatre inconnues : ρ (densité), u (vitesse du fluide), p (pression) et E (énergie totale). L'énergie totale est donnée par :

${\displaystyle E=\rho e+1/2\rho u^{2}}$

où e représente l'énergie interne spécifique.

Pour fermer le système, une équation d'état est nécessaire. Une expression qui convient à notre objectif est celle du gaz parfait :

${\displaystyle p=\rho \left(\gamma -1\right)e}$

où ${\displaystyle \textstyle \gamma }$ est égal au rapport des chaleurs spécifiques ${\displaystyle \textstyle c_{p}/c_{v}}$ du fluide.

Nous pouvons maintenant procéder, comme illustré précédemment dans l'exemple unidimensionnel plus simple, en obtenant les états extrapolés à gauche et à droite pour chaque variable d'état. Ainsi, pour la densité, on obtient :

${\displaystyle \rho _{i+1/2}^{*}=\rho _{i+1/2}^{*}\left(\rho _{i+1/2}^{L},\rho _{i+1/2}^{R}\right),\quad \rho _{i-1/2}^{*}=\rho _{i-1/2}^{*}\left(\rho _{i-1/2}^{L},\rho _{i-1/2}^{R}\right)}$

où

${\displaystyle \rho _{i+1/2}^{L}=\rho _{i}+1/2\,\phi \left(r_{i}\right)\left(\rho _{i+1}-\rho _{i}\right),\quad \rho _{i+1/2}^{R}=\rho _{i+1}-1/2\,\phi \left(r_{i+1}\right)\left(\rho _{i+2}-\rho _{i+1}\right),}$

${\displaystyle \rho _{i-1/2}^{L}=\rho _{i-1}+1/2\,\phi \left(r_{i-1}\right)\left(\rho _{i}-\rho _{i-1}\right),\quad \rho _{i-1/2}^{R}=\rho _{i}-1/2\,\phi \left(r_{i}\right)\left(\rho _{i+1}-\rho _{i}\right)}$

De même, pour la quantité de mouvement ρu et l'énergie totale E, la vitesse u est calculée à partir de la quantité de mouvement et la pression p à partir de l'équation d'état.

Ayant obtenu les états extrapolés limités, nous procédons ensuite à la construction des flux de bord à partir de ces valeurs. Connaissant les flux de bord, nous pouvons alors construire le schéma semi-discret, c'est-à-dire :

${\displaystyle d\mathbf {U} _{i}/dt=-1/\Delta x_{i}\left(\mathbf {F} _{i+1/2}^{*}-\mathbf {F} _{i-1/2}^{*}\right)}$

La solution peut maintenant être obtenue par intégration à l'aide de techniques numériques classiques.

Ce qui précède illustre le principe de base du schéma MUSCL. Cependant, pour une solution pratique aux équations d'Euler, il est également nécessaire de choisir un schéma approprié (tel que le schéma KT présenté ci-dessus) afin de définir la fonction ${\displaystyle \mathbf {F} _{i\pm 1/2}^{*}}$.

Le diagramme ci-contre illustre une solution du second ordre au problème du tube à choc de Sod (en)^[8], obtenue à l'aide du schéma central de Kurganov et Tadmor (KT) haute résolution décrit ci-dessus, avec extrapolation linéaire et limiteur OSPRE. Ceci démontre clairement l'efficacité de l'approche MUSCL pour la résolution des équations d'Euler. La simulation a été réalisée sur un maillage de 200 cellules à l'aide du code MATLAB^[9] adapté pour utiliser l'algorithme KT et le limiteur de flux. L'intégration temporelle a été effectuée par un intégrateur SHK d'ordre 4 (performances équivalentes à celles de RK4). Les conditions initiales suivantes (en unités SI) ont été utilisées :

• pression à gauche = 100 000 Pa
• Pression à droite = 10 000 Pa ;
• Densité à gauche = 1,0 kg/m³ ;
• Densité à droite = 0,125 kg/m³ ;
• Longueur = 20 m ;
• Vitesse à gauche = 0 m/s ;
• Vitesse à droite = 0 m/s ;
• Durée = 0,01 s ;
• λ = 0,001069 (Δt/Δx).

Le diagramme ci-contre illustre une solution du 3e ordre au problème du tube à choc de Sod obtenue à l'aide du schéma central de Kurganov et Tadmor (KT) haute résolution décrit précédemment, mais avec reconstruction parabolique et limiteur de van Albada. Ceci démontre une fois de plus l'efficacité de l'approche MUSCL pour la résolution des équations d'Euler. La simulation a été réalisée sur un maillage de 200 cellules à l'aide dun code MATLAB^[9] adapté pour utiliser l'algorithme KT avec extrapolation parabolique et limiteur de flux (ou limiteur de van Albada). La forme alternative du limiteur de van Albada ${\displaystyle \textstyle \phi _{va}(r)=2r/(1-r^{2})}$ a été utilisée pour éviter les oscillations parasites. L'intégration temporelle a été réalisée à l'aide d'un intégrateur SHK d'ordre 4. Les mêmes conditions initiales ont été utilisées.

Divers autres schémas à haute résolution ont été développés pour résoudre les équations d'Euler avec une bonne précision^[10]. On peut citer, par exemple le schéma de Osher et le schéma AUSM}^[11].