Reading Time: 12 minutes

Lorsque les chercheurs passent de simulations de flux incompressibles à compressibles en Python, ils rencontrent une physique fondamentalement différente – la densité devient une variable, l’équation d’énergie doit être résolue et les ondes de choc introduisent des discontinuités que les méthodes numériques standard ne peuvent tout simplement pas gérer. Cet article couvre ce qui rend SuperSonic Flow différent, les méthodes de capture des chocs conçues pour le gérer et les solveurs Python qui implémentent ces méthodes pour les applications aérospatiales.

Points à retenir clés

  • Le flux compressible nécessite de résoudre les équations compressibles de Navier-Stokes avec la densité comme variable et l’équation énergétique couplée – contrairement au flux incompressible où la densité est constante.
  • Méthodes de capture des chocs (ROE, AUSM, TVD, WENO, Teno) sont essentielles car les chocs introduisent des discontinuités. Le schéma ROE a besoin d’une solution d’entropie et d’une protection contre les carbuncles sur les chocs puissants.
  • Les solveurs de Python pour le flux compressible sont désormais de niveau de production : pybaram (grilles non structurées, 2022), opensbli v3.0 (génération de code automatisée, CPC 2025), JANC (Differentiable GPU Solver, avril 2025) et PygasFlow (quasi-1D).
  • Les cas de repère comme le tube de choc SOD et le problème Shu-Osher fournissent une validation standard pour la précision de la capture des chocs.
  • Muscl-TVD est moins cher mais oscille près des discontinuités ; Weno du 5e ordre donne une résolution supérieure mais coûte plus cher. Pour les DNS/ILES des interactions entre les couches-limites-couches, WENO/TENO est la norme actuelle.
  • Pygasflow et le référentiel open source python-shock-tube fournissent des points d’entrée accessibles pour comprendre la physique des chocs avant de s’engager dans les solveurs de production.

Ce qui différencie le flux supersonique du flux incompressible

En flux incompressible, la densité est traitée comme une constante. Les équations de momentum sont résolues indépendamment de l’équation d’énergie, et la pression est déterminée uniquement par le champ de vitesse par un algorithme de couplage pression-vitesse (comme simple ou piso). Cela signifie que pour les débits à faible vitesse, vous pouvez résoudre les équations de mouvement et de continuité sans jamais connaître la densité du fluide – c’est une donnée.

Dans Compressible Flow, tout change. La densité devient une variable qui évolue en fonction de l’équation de continuité, des équations de mouvement et de l’équation d’énergie. L’état thermodynamique du fluide – pression, température et densité – est couplé par l’équation de l’état (généralement la loi sur les gaz parfaits pour les applications aérospatiales). Cela signifie que vous ne pouvez pas découpler la quantité de mouvement de l’équation de l’énergie. Le solveur doit suivre la propagation des ondes de pression à travers le fluide, et ces ondes se déplacent à la vitesse du son.

Le paramètre sans dimension qui régit cette transition est le nombre de Mach – le rapport de la vitesse d’écoulement à la vitesse locale du son. Quand m < 0,3, les changements de densité sont généralement suffisamment faibles pour que l’hypothèse incompressible le soit. Quand m > 1 (superonique), le flux contient des ondes de choc – des discontinuités où la densité, la pression et la température sautent presque instantanément à travers une région étroite.

Pour les chercheurs en aérospatiale utilisant Python, cela signifie l’amorce Python CFD (post 390) – qui couvre les fluides FIPY, SU2 et JAX principalement pour les flux incompressibles – laisse un écart important. Vous devez comprendre :

  1. Les équations compressibles de Navier-Stokes dans leur forme complète, y compris l’équation de l’énergie
  2. Riemann Solvers qui calcule les flux sur les faces des cellules là où des discontinuités peuvent exister
  3. Méthodes de reconstruction qui déterminent la façon dont les variables s’interpolent des centres cellulaires aux visages
  4. Schémas d’intégration temporelle qui respectent les contraintes de LFC plus strictes à des nombres mach élevés
  5. Traitements numériques spéciaux comme la solution d’entropie et la protection des carbuncles

Il ne s’agit pas d’une simple extension du flux incompressible – c’est une méthodologie numérique différente appliquée à un régime physique différent.

Méthodes de capture des chocs : de ROE à Weno

Lorsqu’une onde de choc se forme, les variables de débit changent de manière discontinue sur une distance inférieure à celle que n’importe quelle grille physique peut résoudre. Au lieu d’essayer de résoudre directement les chocs (ce qui nécessiterait des grilles incroyablement fines), les méthodes de capture des chocs permettent au schéma numérique d’absorber la discontinuité dans quelques cellules, en utilisant la dissipation numérique pour stabiliser la solution.

Le schéma ROE et pourquoi il a besoin de corrections

Le schéma ROE, introduit par Peter Roe en 1978, est un solveur de Riemann approximatif qui linéarise les équations d’Euler autour de l’état de la cellule-face. Il est efficace en termes de calcul et a été la méthode de capture de choc dominante pendant des décennies. L’idée est simple : à chaque face de cellule, vous résolvez une version linéarisée du problème de Riemann (le problème de la valeur initiale pour les données discontinues), calculez le flux et ajoutez une dissipation numérique proportionnelle aux vitesses d’onde.

En Python, la mise en œuvre du schéma ROE à partir de ARXIV 2103.02794 suit ce modèle :

# Roe scheme flux computation (simplified from arXiv 2103.02794)
def roe_flux(uL, uR, n):
    """Compute Roe flux for Euler equations.
    uL, uR: left and right primitive states [rho, u, P]
    n: normal direction (1 or -1 for 1D)
    """
    rhoL, uL, P_L = uL
    rhoR, uR, P_R = uR
    
    # Compute interface states
    rho_hat = (rhoL**2 * rhoR**2 / 
               (rhoL + rhoR))**0.5
    u_hat = (rhoL * uL + rhoR * uR) / (rhoL + rhoR)
    
    # Speed of sound
    S_hat = ((P_L + P_R) / 
             (2 * (rhoL + rhoR)))**0.5
    
    # Compute Roe matrix eigenvalues
    # |u_hat| +/- S_hat for acoustic waves
    # |u_hat| for shear/convection
    # ... eigenvalue computation omitted for brevity

Le problème avec le schéma ROE de base est double. Premièrement, lors de chocs violents, le schéma peut développer une pression négative – un état non physique qui provoque le crash du solveur. Ceci est abordé par la Fix d’entropie, qui modifie les vitesses d’onde à proximité des chocs pour restaurer la convexité de la fonction de flux et prévenir les états non physiques.

Deuxièmement, sur des grilles non structurées, le schéma ROE peut présenter le phénomène du carbuncle – une instabilité numérique où un fort choc (comme celui formé derrière un corps émoussé) développe un deuxième choc non physique ramifié sur le choc principal. Cela se produit parce que le schéma ROE est trop sensible à l’orientation de la grille lors de chocs importants. Le « schéma de ROE stable aux chocs » combine Entropy Fix avec un solveur Riemann tourné pour résoudre ce problème. Pour les chercheurs travaillant sur des géométries aérospatiales avec de forts chocs (comme des cylindres supersoniques ou des corps émoussés), l’omission d’une protection contre les carbonules peut entraîner des défaillances du solveur.

Ausm : fractionner le flux

La méthode d’advection en amont (AUSM) adopte une approche différente. Au lieu de calculer la matrice de ROE, l’AUSM divise le flux en composants d’advection et de pression. Cela le rend robuste à la fois à des nombres mach élevés et à des chocs puissants sans nécessiter de correction d’entropie – la séparation gère intrinsèquement la condition d’entropie physique. Les variantes AUSM (AUSM+, AUSM+, AUSM+) et leurs successeurs (AUSMPWR+) sont parmi les schémas de capture de choc les plus utilisés dans les codes CFD aérospatiaux car ils gèrent la gamme complète de Mach (de Subsonic à hypersonique) sans changement de stratégie.

Le référentiel Python-Shock-Tube implémente AUSM avec une interface propre qui fonctionne aux côtés de Godounov, Roe, HLL et HLLC :

# AUSM implementation pattern from Python-shock-tube
def ausm_flux(uL, uR, n):
    """AUSM flux with pressure splitting."""
    rhoL, uL, P_L = uL
    rhoR, uR, P_R = uR
    
    # Split flux into advective and pressure parts
    # Advective part uses Mach number-based splitting
    # Pressure part uses sign function of Mach number
    # See Toro's "Riemann Solvers and Numerical Methods 
    # for Fluid Dynamics" for full derivation

Limiteurs de flux TVD et reconstruction MUSCL

Les limiteurs de flux TVD (diminution totale de la variation) empêchent les oscillations parasites à proximité des discontinuités en contrôlant la variation totale de la solution. Lorsqu’ils sont combinés à la reconstruction de MUSCL (monotonic Upstream-Centered for Conservation Laws), ils forment le schéma MusCl-TVD – une méthode de second ordre qui est la bourreau de travail pour la capture de choc au niveau technique.

Le compromis ici est clair : MUSCL-TVD est peu coûteux et robuste, mais à des chocs violents ou à des interactions avec des ondes complexes, il peut produire de petites oscillations derrière le front de choc. Pour les applications d’ingénierie où la résolution absolue des chocs n’est pas critique, MUSCL-TVD est souvent suffisant.

Weno et Teno : capture de chocs d’ordre élevé

La reconstruction pondérée essentiellement non oscillatoire (Weno) du 5e ordre aborde les limites de MUSCL-TVD en utilisant une combinaison pondérée de pochoirs qui s’adaptent automatiquement aux régions lisses et discontinues. Près d’un débit doux, tous les pochoirs contribuent ; Près des chocs, les poids se concentrent sur le pochoir non oscillatoire. Le résultat est une résolution considérablement meilleure des oscillations à haute fréquence proches des discontinuités – cruciales pour la simulation numérique directe haute fidélité (DNS) et la simulation implicite des grands Feddy (ILES).

Le référentiel Python-Shock-Tube le démontre clairement avec le Shu–Osher Problem – une référence classique où un choc interagit avec une onde de densité. Voir l’image de comparaison Shu–Osher du référentiel, qui montre que WENO du 5e ordre résout bien mieux les oscillations à haute fréquence derrière le choc que MUSCL-TVD.

Pour les DNS et les iles d’interactions entre les couches limite de choc – comme les cas de buffet transonique validés dans OpenSBLI v3.0 – Weno et sa variante teno préservant les bords (enO ciblé) sont la norme actuelle. La génération de code automatisée d’OpenSBLI produit du code C de qualité de production à partir de formulations symboliques Python utilisant Weno/Teno pour la discrétisation compressible de Navier-Stokes.

Solveurs Python pour un flux aérospatial compressible

Le paysage des solveurs de flux compressibles basés sur Python a considérablement mûri entre 2021 et 2025. Cinq cadres sont désormais disponibles pour les chercheurs en aérospatiale, chacun avec des architectures distinctes :

Pybaram — Grilles non structurées pour les chocs forts

Pybaram (publié dans SoftwareX, 2022) met en œuvre des méthodes de volume fini sur des grilles non structurées pour les géométries aérospatiales 2D et 3D. Il gère explicitement les fortes ondes de choc et d’expansion, conçues pour l’enseignement et la recherche avec des données ouvertes sur Mendeley. Avec plus de 10 citations, Pybaram représente un point d’entrée accessible pour les chercheurs qui ont besoin d’une prise en charge non structurée de la grille sans les frais généraux de plus grands packages.

OpenSBLI v3.0 — génération de code automatisée pour DNS/ILES

OpenSbli représente une approche fondamentalement différente : au lieu d’écrire directement du code CFD compilé, les utilisateurs écrivent des équations dans la notation Einstein à l’aide de l’algèbre symbolique Python (Sympy) et le générateur de code d’OpenSbli produit un code C optimisé via l’OPS DSL. Le framework cible les grilles structurées multi-blocs et prend en charge les backends CPU (MPI/OpenMP) et GPU (CUDA/OpenCL).

OpenSbli v3.0, publié dans Ordinateurs & Physics (janvier 2025), prend en charge explicitement les Weno et Teno Shock Capture et a été validé sur :

  • Interactions de la couche limite de choc au niveau des nombres de Mach hypersoniques
  • Buffet transonique sur les profils aérodynamiques de la NASA-CRM
  • DNS du flux hypersonique sur les bosses gaussiennes

Le pedigree académique du cadre est solide – l’équipe de Southampton/Oxford (Lusher, Sandham) a publié plus de 20 articles sur les applications OpenSBLI. Si vos recherches portent sur des DN/ILES haute fidélité de flux compressibles avec interactions de choc, la génération de code automatisée d’OpenSBLI est l’option la plus mature.

Le diagramme de flux de travail d’OpenSBL Framework illustre comment les équations symboliques de Python sont traduites en code C optimisé via l’OPS DSL.

JANC — Flux compressible différentiable sur GPU

JANC (publié en avril 2025 par des chercheurs de l’Université de Tsinghua et des collaborateurs) est sans doute le solveur compressible Python le plus novateur de ces dernières années. Construit entièrement sur JAX-AMR (Raffinement de maillage adaptatif structuré basé sur JAX), JANC est un solveur entièrement différentiable pour les flux réactifs compressibles. Sa prétention à la distinction :

Sur un seul GPU A100, JANC réalise 1 % du coût de calcul d’Openfoam (utilisant 384 cœurs). Avec l’AMR activée, le coût chute à 1 à 2 % d’OpenFoam.

Il ne s’agit pas d’une allégation marketing – elle est vérifiée par rapport aux cas de test signalés. Leviers de JANC :

  • Différenciation automatique (pas besoin de dériver manuellement les équations adjointes)
  • Compilation XLA JIT (optimisation du calcul du tenseur de XLA)
  • Raffinement de maillage adaptatif (résolution automatique où les gradients sont raides)
  • Parallélisme natif GPU/TPU
  • Licence MIT

Pour les chercheurs qui explorent les CFD différenciables – une tendance croissante pour l’optimisation de la conception, l’analyse de sensibilité et l’intégration de ML – JANC est le premier solveur python entièrement différenciable pour réagir des flux compressibles. Le code est disponible sous licence MIT.

PygasFlow et SU2 — Calculs rapides et suites aérospatiales complètes

PygasFlow (v1.4.1) est un ensemble léger Python pour la dynamique des gaz parfaits quasi-1D. Il met en œuvre des relations de choc normales et obliques, un flux isentropique, un flux de Fanno et Rayleigh et des corrélations aérothermodynamiques de plaques plates/stagnantes. Bien qu’il ne s’agisse pas d’un solveur CFD complet, il est inestimable pour des calculs analytiques rapides – diagrammes interactifs de déflexion de pression, contrôles de relation entre les chocs et balayages paramétriques rapides avant de s’engager dans une simulation complète.

SU2 reste la suite CFD open-source la plus citée pour les applications aérospatiales. Ses capacités supersoniques/compressibles comprennent le Wedge Supersonic Supersonic, le profil aérodynamique Transonic OneRam6 et le flux compressible non idéal dans les didacticiels de buses Supersonic. Avec plus de 4 400 étoiles GitHub, SU2 est le choix établi pour la production aérospatiale CFD où Python sert principalement de couche de script et d’automatisation.

Cas de référence et validation

Toute implémentation de capture de chocs nécessite une validation par rapport aux cas de référence standard. Deux cas sont essentiels :

Tube de choc de gazon

Le tube de choc SOD est le test canonique pour les solveurs Riemann. Il initialise un tube avec un état gauche haute pression et un état de droite basse pression. À t = 0, le diaphragme se brise et un choc, une discontinuité de contact et une forme de ventilateur d’expansion. Le référentiel python-shock-tube fournit cette référence. Voir la comparaison de tube de choc SOD montrant la résolution du flux sur quatre solveurs (Godunov, ROE, HLL, HLLC) sur 300 cellules avec reconstruction d’ordre 0.

L’image montre une résolution de flux sur quatre solveurs (Godunov, ROE, HLL, HLLC) sur 300 cellules avec une reconstruction d’ordre 0. HLLC (hybride Linéarisé LAX – Friedrich avec contact) fournit généralement le meilleur équilibre : il résout mieux la discontinuité de contact que HLL tout en étant comparable à HLL. C’est pourquoi HLLC est le choix par défaut dans l’implémentation WENO du référentiel python-shock-tube.

Shu-Osher Problème

Le problème de Shu-Osher teste dans quelle mesure un schéma gère l’interaction entre un choc et une onde de densité. La condition initiale est une fonction d’étape pour la pression avec une perturbation de densité sinusoïdale derrière le choc. Ce cas révèle comment les méthodes de reconstruction résolvent les oscillations à haute fréquence :

  • La reconstruction du 0e ordre (Godunov) subit une dissipation numérique excessive lors de chocs violents
  • Muscl-TVD réduit la dissipation mais oscille derrière le choc
  • Weno du 5e ordre résout avec précision les oscillations à haute fréquence

Pour les simulations aérospatiales haute fidélité, Weno est désormais le choix standard.

Défis numériques uniques à un flux mach élevé

Plusieurs défis numériques surviennent spécifiquement dans le flux compressible/supsonique qui n’apparaissent pas (ou n’apparaissent pas différemment) dans le CFD incompressible :

Contraintes de LFC à des nombres mach élevés

La condition de LCF (Courant – Friedrichs-Lewy) régit la stabilité pour l’intégration explicite dans le temps. Pour un écoulement compressible, le nombre de LCF dépend de la somme de la vitesse d’écoulement et de la vitesse du son :

$$CFL = frac{u + a}{delta x} delta t$$

À des nombres mach élevés, $u approx m cdot a$, de sorte que la condition LFC devient extrêmement restrictive. Cela signifie que pas de temps sont requis, souvent des ordres de grandeur inférieurs à ceux des flux incompressibles. Les méthodes implicites de pas de temps ou de multi-grilles peuvent atténuer cela, mais les schémas explicites nécessitent un contrôle attentif du pas de temps.

Correction d’entropie

Le schéma ROE de base échoue à la condition d’entropie lors de chocs violents : au lieu de l’augmentation d’entropie physiquement correcte à travers le choc, la solution numérique peut violer la thermodynamique. La correction d’entropie modifie les valeurs propres à proximité des chocs pour garantir que le schéma numérique respecte la deuxième loi de la thermodynamique. Sans cela, la solution peut converger vers un mauvais état d’entropie – un choc trop fort ou trop faible.

phénomène de l’escarboulement

Comme mentionné précédemment, le phénomène de carbuncle est une instabilité numérique spécifique au schéma ROE sur des grilles non structurées. Lorsqu’un fort choc s’aligne mal avec la grille, le solveur peut développer un choc secondaire qui se démarque du choc principal – physiquement impossible et numériquement catastrophique. Le schéma de ROE stable aux chocs (combinant le correctif d’entropie avec le solveur de Riemann pivoté) guérit cela.

Interaction de la couche limite choc-limite

Dans les applications aérospatiales, les ondes de choc interagissant avec les couches limites produisent des structures d’écoulement complexes – séparation des flux, rattachement et buffet instable. La résolution de ces éléments nécessite :

  • Résolution fine de la couche limite (souvent DNS résolu en mur)
  • Capture de choc d’ordre élevé (Weno/Teno) pour résoudre le choc
  • Résolution temporelle adéquate pour les phénomènes instables

C’est pourquoi la validation d’OpenSbli v3.0 se concentre sur SBLI – c’est la physique qui compte le plus pour la conception de véhicules transoniques et hypersoniques.

Quel solveur choisir ?

Le choix du solveur de flux compressible dépend des besoins de votre réseau, des besoins de fidélité et des ressources de calcul. Voici une comparaison :

solveur Type de grille capture des chocs le mieux pour Matériel
pybaram non structuré volume fini Enseignement, 2D/3D aérospatial, chocs puissants CPU parallèle
opensbli Structuré (multibloc) Weno/Teno DNS/ILES, SBLI, buffet transonique, hypersonique CPU/GPU (CUDA/OpenCL)
janc Structuré (Bloc AMR) Jax Auto-Diff CFD différentiable, optimisation adjointe, GPU GPU/TPU (A100+)
Pygasflow N/A (analytique) N/A Relations de choc rapide, Fanno/Rayleigh, balayages paramétriques Tout
su2 Structuré/non structuré ROE, Ausm, etc. Suite aérospatiale complète, CFD de production, flux N+1 CPU/MPI

Recommandation pratique : commencez par PygasFlow pour la compréhension analytique des relations de choc et du flux de Fanno/Rayleigh. Utilisez ensuite le référentiel Python-Shock-Tube pour la pratique de la capture de chocs 1D (ROE, AUSM, WENO). Choisissez Pybaram si vous avez besoin de grilles non structurées pour les géométries complexes 2D/3D. Choisissez OpenSBLi si votre recherche implique DNS/ILES avec WENO/TENO sur des grilles structurées. Choisissez JANC si vous explorez CFD différentiable avec GPU Acceleration et AMR.

Résumé et étapes suivantes

Cet article couvrait ce qui rend le flux compressible/superonique différent du flux incompressible, les méthodes de capture des chocs (ROE avec entropy Fix, AUSM, TVD, WENO/TENO) qui gèrent les discontinuités, les solveurs Python disponibles pour les applications aérospatiales (Pybaram, OpenSBLI, JANC, PygasFlow, SU2) et cas de référence pour la validation.

Pour les chercheurs commençant par un flux compressible en Python : commencez par PygasFlow pour des calculs analytiques rapides, pratiquez des méthodes de capture de chocs avec Python-Shock-Tube, puis choisissez Pybaram (non structuré) ou OpenSbli (structuré) en fonction des exigences de votre grille. Pour les CFD différenciables avec accélération GPU, JANC représente l’état actuel de la technique.

Si vous débutez dans Python CFD, commencez par CFD complet couvrant Fipy, maillage et les 12 étapes de Navier–Stokes. Pour les flux dominés par l’advection où la résolution de maillage est importante, Adaptation du maillage anisotrope peut compléter votre schéma de capture de chocs. Lorsque vous avez besoin d’une meilleure résolution que le volume fini de second ordre, les méthodes spectrales et DG d’ordre élevé offrent des alternatives.

Lecture recommandée