Prérequis : Chapitre 3 — Méthode des volumes finis (semi-discrétisation, cellules fantômes, principe de composition reconstruction + flux), Chapitre 4 — Flux numériques (solveurs de Riemann composables).
Un schéma de volumes finis d'ordre 1 utilise une reconstruction constante par morceaux : dans chaque cellule
Peut-on simplement utiliser une reconstruction linéaire (ou polynomiale de degré supérieur) sans précaution particulière ? Non. Le théorème de Godunov (1959) affirme :
Un schéma linéaire monotone (ne créant pas de nouveaux extrema) est au plus d'ordre 1.
Ce résultat fondamental signifie que toute montée en ordre passe nécessairement par une reconstruction non linéaire, c'est-à-dire dont les coefficients dépendent de la solution locale. C'est exactement ce que font les méthodes MUSCL (limiteurs de pente), ENO (sélection de stencil) et WENO (pondération non linéaire de stencils) présentées dans ce chapitre.
Toutes les reconstructions de ce solveur opèrent sur les variables primitives
- Près d'un choc,
$\rho u$ et$E$ présentent des sauts couplés (le produit$\rho u$ mélange densité et vitesse), tandis que$\rho$ ,$u$ et$p$ ont des sauts indépendants et plus réguliers. - Près d'une discontinuité de contact, seule
$\rho$ est discontinue parmi les primitives, alors que$\rho u$ et$E$ le sont aussi parmi les conservatives. - Les limiteurs et indicateurs de régularité fonctionnent mieux sur des variables dont les variations sont découplées et monotones par morceaux.
Le flux de calcul est donc :
Un schéma d'ordre élevé se compose de deux briques indépendantes (cf. §3.5) :
La reconstruction fournit les états
Toutes les reconstructions de ce solveur appliquent un plancher de positivité aux valeurs reconstruites de
Ce garde-fou empêche l'apparition de densités ou pressions négatives qui provoqueraient des erreurs numériques (racine carrée d'un nombre négatif dans le calcul de la vitesse du son).
La méthode MUSCL (Monotone Upstream-centered Scheme for Conservation Laws), introduite par van Leer [van Leer, 1979], est la plus répandue des reconstructions d'ordre 2. L'idée est simple : au lieu de représenter la solution par une constante dans chaque cellule, on utilise une droite (reconstruction affine). La pente de cette droite est limitée pour respecter une condition de type TVD (Total Variation Diminishing), ce qui empêche la création de nouveaux extrema.
Dans chaque cellule
où
On définit les différences forward aux interfaces :
Le rapport de pente dans la cellule
La pente limitée est alors :
où
Pour qu'un limiteur
-
$\varphi(r) = 0$ pour$r \leq 0$ (annulation de la pente aux extrema). - La courbe
$\varphi(r)$ doit se situer dans la région de Sweby, délimitée par :-
Borne inférieure : le limiteur minmod,
$\varphi(r) = \min(1, r)$ . -
Borne supérieure : le limiteur superbee,
$\varphi(r) = \max(\min(2r, 1), \min(r, 2))$ .
-
Borne inférieure : le limiteur minmod,
-
$\varphi(1) = 1$ (reconstruction centrée quand la solution est localement linéaire).
Tout limiteur dont le graphe reste dans cette région est TVD et d'ordre 2 en zone lisse.
- Propriétés : limiteur le plus dissipatif de la région TVD (borne inférieure de Sweby). Sélectionne la pente de plus petite valeur absolue entre la différence amont et la différence aval.
- Avantages : très robuste, jamais d'oscillation, convergence garantie. Idéal pour les problèmes raides.
- Inconvénients : trop dissipatif en zone lisse ; lisse excessivement les discontinuités de contact et les gradients doux. Réduit l'ordre effectif près des extrema lisses.
-
Propriétés : limiteur différentiable (lisse en
$r = 0$ ), situé au milieu de la région de Sweby. Bon compromis entre dissipation et précision. - Avantages : pas de discontinuité dans la fonction limiteur, ce qui favorise la convergence des solveurs implicites. Excellent compromis dissipation/précision pour la plupart des applications.
- Inconvénients : légèrement plus dissipatif que superbee ou MC sur les contacts.
- Propriétés : limiteur le moins dissipatif de la région TVD (borne supérieure de Sweby). Maximise la pente dans les limites TVD.
- Avantages : excellent pour les discontinuités de contact qu'il maintient très raides. Très faible dissipation numérique.
- Inconvénients : peut artificiellement raidir des profils lisses (effet de « crénelage » ou staircase effect), transformant des transitions douces en marches d'escalier. Ce comportement compressif est parfois indésirable.
-
Propriétés : utilise la pente centrée
$(1+r)/2$ contrainte par les bornes TVD$2r$ et$2$ . Se situe entre van Leer et superbee. - Avantages : meilleure précision que van Leer en zone lisse grâce à la pente centrée. Bonne résolution des contacts sans l'effet compressif de superbee.
- Inconvénients : non différentiable (comme minmod et superbee), ce qui peut ralentir la convergence des schémas implicites.
-
Propriétés : limiteur différentiable utilisant un mélange quadratique rationnel. Comme van Leer, la transition est lisse en
$r = 0$ et$r = 1$ . Se situe dans la partie basse de la région de Sweby, entre minmod et van Leer. -
Avantages : continuité
$C^1$ qui favorise la convergence des schémas implicites. Produit des solutions très régulières. Satisfait la propriété TVD. -
Inconvénients : plus dissipatif que van Leer, surtout près des discontinuités de contact. Les termes quadratiques
$r^2$ au numérateur et dénominateur atténuent davantage les gradients modérés.
Référence : [van Albada et al., 1982] G. D. van Albada, B. van Leer, W. W. Roberts, "A comparative study of computational methods in cosmic gas dynamics", Astronomy & Astrophysics, 108, pp. 76–84.
| Limiteur | Formule | Dissipation | Différentiable | Position Sweby | Risque compressif |
|---|---|---|---|---|---|
| Minmod | Maximale | Non | Borne inférieure | Aucun | |
| Van Albada | Élevée | Oui | Bas (entre minmod et van Leer) | Aucun | |
| Van Leer | Modérée | Oui | Centre | Aucun | |
| MC | Faible | Non | Centre-haut | Faible | |
| Superbee | Minimale | Non | Borne supérieure | Élevé |
- Ordre spatial : 2.
-
Cellules fantômes :
$n_{\text{ghost}} = 2$ (une différence de chaque côté du stencil MUSCL). - Intégrateur temporel par défaut : RK2 (Heun), cohérent avec l'ordre spatial 2.
- [van Leer, 1979] B. van Leer, Towards the Ultimate Conservative Difference Scheme. V. A Second-Order Sequel to Godunov's Method, J. Comput. Phys., 32, pp. 101–136.
- [Sweby, 1984] P. K. Sweby, High Resolution Schemes Using Flux Limiters, SIAM J. Numer. Anal., 21(5), pp. 995–1011.
La méthode ENO (Essentially Non-Oscillatory), introduite par Harten, Engquist, Osher et Chakravarthy [Harten et al., 1987], propose une approche radicalement différente des limiteurs de pente : au lieu de limiter une pente fixe, on choisit le stencil le plus lisse parmi plusieurs candidats. L'idée est d'éviter d'interpoler à travers une discontinuité en sélectionnant automatiquement le stencil qui ne la traverse pas.
Pour une reconstruction d'ordre 2 (ENO2), on dispose de deux stencils candidats pour calculer la pente dans la cellule
-
Stencil gauche :
${i-1, i}$ utilisant la différence$D_{i-1} = W_i - W_{i-1}$ . -
Stencil droit :
${i, i+1}$ utilisant la différence$D_i = W_{i+1} - W_i$ .
Les premières différences (ou différences finies forward) sont :
Les secondes différences (différences divisées d'ordre 2) sont :
Le stencil est choisi en comparant les secondes différences de part et d'autre :
L'idée est que la seconde différence la plus petite en valeur absolue correspond à la zone la plus lisse.
Pour l'état gauche à l'interface
où
Pour l'état droit à l'interface
où
Avantages :
- Pas de paramètre à régler (contrairement aux limiteurs qui offrent un choix entre plusieurs fonctions
$\varphi$ ). - Sélection automatique du stencil basée sur un critère objectif (régularité locale).
- Bonne capture des discontinuités sans oscillation.
Inconvénients :
- La commutation abrupte entre stencils peut introduire de petites perturbations quand la solution est proche du seuil de décision.
- Les poids ne sont pas des fonctions lisses de la solution (contrairement à WENO).
- Limité à l'ordre 2 dans cette implémentation.
- Ordre spatial : 2.
-
Cellules fantômes :
$n_{\text{ghost}} = 2$ . - Intégrateur temporel par défaut : RK2 (Heun).
- [Harten et al., 1987] A. Harten, B. Engquist, S. Osher, S. R. Chakravarthy, Uniformly High Order Accurate Essentially Non-Oscillatory Schemes, III, J. Comput. Phys., 71, pp. 231–303.
La méthode ENO sélectionne un seul stencil et rejette les autres, ce qui constitue un gaspillage d'information. De plus, la commutation binaire entre stencils est non lisse. La méthode WENO (Weighted ENO), introduite par Liu, Osher et Chan (1994) puis formalisée par Jiang et Shu [Jiang, Shu, 1996], corrige ces deux défauts en combinant tous les stencils avec des poids non linéaires qui s'adaptent à la régularité locale.
En zone lisse, les poids tendent vers leurs valeurs optimales (qui maximisent l'ordre de précision). Près d'une discontinuité, le poids du stencil qui traverse la discontinuité tend vers zéro, reproduisant le comportement ENO.
Pour WENO3, on utilise deux sous-stencils de 2 cellules chacun :
-
$S_0 = {i-1, i}$ : stencil biaisé à gauche. -
$S_1 = {i, i+1}$ : stencil biaisé à droite.
Chaque stencil fournit une reconstruction linéaire de la solution.
Les indicateurs de régularité (smoothness indicators) mesurent la variation locale sur chaque sous-stencil. Pour WENO3, ils sont simplement le carré des premières différences :
Un indicateur
Les poids idéaux (ou linear weights) sont ceux qui, en zone parfaitement lisse, donnent l'ordre de précision maximal. Pour la reconstruction de l'état gauche :
Pour la reconstruction de l'état droit :
Les poids non linéaires de Jiang-Shu [Jiang, Shu, 1996] sont :
avec
Comportement :
- En zone lisse (
$\beta_0 \approx \beta_1$ ), les poids convergent vers les poids idéaux$d_k$ . - Près d'une discontinuité (
$\beta_k \gg \beta_l$ ), le poids$\omega_k \to 0$ : le stencil contaminé est exclu.
Limitation : aux points critiques (où les dérivées s'annulent,
Pour corriger la perte de précision aux points critiques, Borges et al. [Borges et al., 2008] proposent un indicateur global :
Les poids WENO-Z sont :
avec
Avantage : aux points critiques,
L'état gauche à l'interface
L'état droit se reconstruit de manière symétrique depuis la cellule
- Ordre spatial : 3 (en zone lisse).
-
Cellules fantômes :
$n_{\text{ghost}} = 2$ . - Intégrateur temporel par défaut : SSP-RK3 (Shu-Osher), préservant la propriété TVD.
- [Jiang, Shu, 1996] G.-S. Jiang, C.-W. Shu, Efficient Implementation of Weighted ENO Schemes, J. Comput. Phys., 126, pp. 202–228.
- [Borges et al., 2008] R. Borges, M. Carmona, B. Costa, W. S. Don, An Improved Weighted Essentially Non-Oscillatory Scheme for Hyperbolic Conservation Laws, J. Comput. Phys., 227, pp. 3191–3211.
WENO5 étend le principe de WENO3 en utilisant trois sous-stencils de 3 cellules chacun, atteignant l'ordre 5 en zone lisse. C'est le schéma WENO le plus utilisé en pratique pour les lois de conservation hyperboliques, offrant un excellent compromis entre précision élevée et robustesse près des discontinuités.
Pour la reconstruction de l'état gauche à l'interface
-
$S_0 = {i, i+1, i+2}$ : stencil biaisé à droite. -
$S_1 = {i-1, i, i+1}$ : stencil central. -
$S_2 = {i-2, i-1, i}$ : stencil biaisé à gauche.
Les indicateurs de régularité pour WENO5 sont plus complexes que pour WENO3, car ils intègrent les termes de dérivée seconde du polynôme d'interpolation :
Le premier terme de chaque
Pour la reconstruction de l'état gauche :
Pour la reconstruction de l'état droit, les poids sont inversés :
Le stencil central
avec
L'indicateur global pour WENO5 est :
Les poids sont :
avec
L'indicateur
Chaque sous-stencil fournit un polynôme d'interpolation de degré 2. Pour l'état gauche à l'interface
Pour l'état droit à l'interface
La valeur reconstruite finale est la combinaison pondérée :
| Propriété | JS | Z |
|---|---|---|
| Précision en zone lisse | Ordre 5 | Ordre 5 |
| Précision aux points critiques ( |
Dégradée (ordre 3-4) | Restaurée (ordre 5) |
| Paramètre |
|
|
| Coût supplémentaire | — | Calcul de |
En pratique, la variante Z est recommandée pour les solutions contenant des extrema lisses (ondes acoustiques, tourbillons) car elle préserve l'ordre de convergence théorique.
- Ordre spatial : 5 (en zone lisse).
-
Cellules fantômes :
$n_{\text{ghost}} = 3$ (stencils de 5 cellules). - Intégrateur temporel par défaut : Dormand-Prince RK5 (6 étages), pour que l'erreur temporelle n'entame pas la précision spatiale d'ordre 5.
- [Jiang, Shu, 1996] G.-S. Jiang, C.-W. Shu, Efficient Implementation of Weighted ENO Schemes, J. Comput. Phys., 126, pp. 202–228.
- [Borges et al., 2008] R. Borges, M. Carmona, B. Costa, W. S. Don, An Improved Weighted Essentially Non-Oscillatory Scheme for Hyperbolic Conservation Laws, J. Comput. Phys., 227, pp. 3191–3211.
| Propriété | MUSCL | ENO2 | WENO3 | WENO5 |
|---|---|---|---|---|
| Ordre spatial | 2 | 2 | 3 | 5 |
|
Cellules fantômes ( |
2 | 2 | 2 | 3 |
| Nombre de stencils | 1 (pente limitée) | 2 (sélection) | 2 (pondération) | 3 (pondération) |
| Contrôle des oscillations | Limiteur TVD | Sélection de stencil | Poids non linéaires | Poids non linéaires |
| Précision aux points critiques (JS) | — | — | Dégradée | Dégradée |
| Précision aux points critiques (Z) | — | — | Restaurée | Restaurée |
| Coût relatif | Faible | Faible | Modéré | Élevé |
| Intégrateur par défaut | RK2 | RK2 | SSP-RK3 | RK5 (Dormand-Prince) |
| Paramètres utilisateur | Choix du limiteur | Aucun | Variante JS/Z | Variante JS/Z |
| Composable avec | Rusanov, HLL, HLLC, Roe, Godunov | idem | idem | idem |
Recommandations pratiques :
- Pour un premier calcul rapide ou un problème dominé par des chocs forts : MUSCL + van Leer offre le meilleur compromis robustesse/coût.
- Pour une étude de convergence en maillage ou une solution lisse : WENO5-Z préserve l'ordre 5 même aux extrema.
- Pour une capture fine des contacts : MUSCL + superbee (avec prudence sur les profils lisses) ou WENO5.
- Le choix du flux numérique (HLLC recommandé) a souvent plus d'impact que le choix de la reconstruction sur la qualité de la discontinuité de contact.