Construire un solveur Navier-Stokes en Python à partir de zéro : simulation du flux d’air

(CFD) est souvent considéré comme une boîte noire de logiciels commerciaux complexes. Cependant, implémenter un solveur « à partir de zéro » constitue l’un des moyens les plus puissants d’apprendre la physique du mouvement des fluides. J’ai commencé cela comme un projet personnel, et dans le cadre d’un cours de biophysique, j’en ai profité pour enfin comprendre comment fonctionnent ces belles simulations.

Ce guide est conçu pour les data scientists et les ingénieurs qui souhaitent aller au-delà des bibliothèques de haut niveau et comprendre les mécanismes sous-jacents des simulations numériques en traduisant des équations aux dérivées partielles en code Python discrétisé. Nous explorerons également des concepts fondamentaux de programmation tels que les opérations vectorisées avec NumPy et la convergence stochastique, qui sont des compétences essentielles pour toute personne intéressée par des architectures plus larges de calcul scientifique et d’apprentissage automatique.

Nous passerons en revue la dérivation et l’implémentation Python d’un solveur Navier-Stokes (NS) simple et incompressible. Ensuite, nous appliquerons ce solveur pour simuler le flux d’air autour du profil de l’aile d’un oiseau.

La physique : Navier-Stokes incompressible

Les équations fondamentales du CFD sont les équations de Navier-Stokes, qui décrivent l’évolution de la vitesse et de la pression dans un fluide. Pour un vol régulier (comme un oiseau planant), nous supposons que l’air est incompressible (densité constante) et laminaire. Ils peuvent être compris comme simplement la loi du mouvement de Newton, mais pour un élément infinitésimal de fluide, avec les forces qui l’affectent. Il s’agit principalement de pression et de viscosité, mais selon le contexte, vous pouvez ajouter de la gravité, des contraintes mécaniques ou même de l’électromagnétisme si vous en avez envie. L’auteur peut attester que ce n’est vraiment pas recommandé pour un premier projet.

Les équations sous forme vectorielle sont :

\[
\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{v} \\
\nabla \cdot \mathbf{v} = 0
\]

Où:

• v: Champ de vitesse (u,v)
• p: Pression
• ρ: Densité du fluide
• ν: Viscosité cinématique

La première équation (Momentum) équilibre l’inertie avec les gradients de pression et la diffusion visqueuse. La deuxième équation (Continuité) impose que la densité du fluide reste constante.

Le problème du couplage de pression

Un défi majeur en CFD est que la pression et la vitesse sont couplées : le champ de pression doit s’ajuster constamment pour garantir que le fluide reste incompressible.

Pour résoudre cela, nous dérivons un Équation de pression-Poisson en prenant la divergence de l’équation de la quantité de mouvement. Dans un solveur discrétisé, nous résolvons cette équation de Poisson à chaque pas de temps pour mettre à jour la pression, garantissant ainsi que le champ de vitesse reste sans divergence.

Discrétisation : des mathématiques à la grille

Pour résoudre ces équations sur un ordinateur, nous utilisons Différence finie schémas sur une grille uniforme.

• Temps: Différence à terme (Euler explicite).
• Advection (termes non linéaires) : Différence arrière/montée au vent (pour la stabilité).
• Diffusion et pression : Différence centrale.

Par exemple, la formule de mise à jour pour le composant u (vitesse x) ressemble à ceci sous forme de différences finies

\[
u_{i,j}^n \frac{u_{i,j}^n - u_{i-1,j}^n}{\Delta x}
\]

Dans le code, le terme d’advection u∂x∂u​ utilise une différence vers l’arrière :

\[
u_{i,j}^n \frac{u_{i,j}^n - u_{i-1,j}^n}{\Delta x}
\]

L’implémentation Python

L’implémentation se déroule en quatre étapes distinctes à l’aide de tableaux NumPy.

1. Initialisation

Nous définissons la taille de la grille (nx, ny), pas de temps (dt) et les paramètres physiques (rho, nu). Nous initialisons les champs de vitesse (u,v) et de pression (p) à des zéros ou à un flux uniforme.

2. La géométrie de l’aile (limite immergée)

Pour simuler une aile sur une grille cartésienne, nous devons marquer quels points de la grille se trouvent à l’intérieur l’aile solide.

• Nous chargeons un maillage d’aile (par exemple, à partir d’un fichier STL).
• Nous créons un tableau de masques booléens où True indique un point à l’intérieur de l’aile.
• Lors de la simulation, nous forçons la vitesse à zéro en ces points masqués (condition de non-glissement/non-pénétration).

3. La boucle principale du solveur

La boucle centrale se répète jusqu’à ce que la solution atteigne un état stable. Les étapes sont les suivantes :

1. Construisez le terme source (b) : Calculez la divergence des termes de vitesse.
2. Résoudre la pression : Résolvez l’équation de Poisson pour p en utilisant l’itération de Jacobi.
3. Vitesse de mise à jour : Utilisez la nouvelle pression pour mettre à jour u et v.
4. Appliquer des conditions aux limites : Appliquez une vitesse d’entrée et des vitesses nulles à l’intérieur de l’aile.

Le code

Voici à quoi ressemblent les mises à jour mathématiques de base en Python (vectorisées pour les performances).

Étape A : Création du terme source de pression Cela représente le côté droit (RHS) de l’équation de Poisson basée sur les vitesses actuelles.

# b is the source term
# u and v are current velocity arrays
b[1:-1, 1:-1] = (rho * (
    1 / dt * ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx) +
              (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
    ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))**2 -
    2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
         (v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx)) -
    ((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))**2
))

Étape B : Résolution de la pression (itération Jacobi) Nous itérons pour lisser le champ de pression jusqu’à ce qu’il équilibre le terme source.

for _ in range(nit):
    pn = p.copy()
    p[1:-1, 1:-1] = (
        (pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +
        (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2 -
        b[1:-1, 1:-1] * dx**2 * dy**2
    ) / (2 * (dx**2 + dy**2))
# Boundary conditions: p=0 at edges (gauge pressure)
    p[:, -1] = 0; p[:, 0] = 0; p[-1, :] = 0; p[0, :] = 0

Étape C : mise à jour de la vitesse Enfin, nous mettons à jour la vitesse en utilisant les équations de moment discrétisées explicites.

un = u.copy()
vn = v.copy() 
# Update u (x-velocity)
u[1:-1, 1:-1] = (un[1:-1, 1:-1] -
                 un[1:-1, 1:-1] * dt / dx * (un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
                 vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
                 dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, 0:-2]) +
                 nu * (dt / dx**2 * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
                       dt / dy**2 * (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])))
# Update v (y-velocity)
v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
                 un[1:-1, 1:-1] * dt / dx * (vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
                 vn[1:-1, 1:-1] * dt / dy * (vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
                 dt / (2 * rho * dy) * (p[2:, 1:-1] - p[0:-2, 1:-1]) +
                 nu * (dt / dx**2 * (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
                       dt / dy**2 * (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))

Résultats : Est-ce que ça vole ?

Nous avons exécuté ce solveur sur un profil d’aile rigide avec un afflux constant en champ lointain.

Observations qualitatives Les résultats correspondent aux attentes physiques. Les simulations montrent une pression élevée sous l’aile et une pression faible au-dessus, ce qui est exactement le mécanisme qui génère la portance. Les vecteurs de vitesse montrent l’accélération du flux d’air sur la surface supérieure (principe de Bernoulli).

Forces : soulèvement ou traînée En intégrant le champ de pression sur la surface de l’aile, nous pouvons calculer la portance.

• Le solveur démontre que les forces de pression dominent les forces de friction visqueuses sont multipliées par près de 1 000 dans l’air.
• À mesure que l’angle d’attaque augmente (de 0∘ à −20∘), le rapport portance/traînée augmente, correspondant aux tendances observées dans les souffleries et les logiciels CFD professionnels comme OpenFOAM.

Limites et prochaines étapes

Bien que la création de ce solveur ait été idéale pour l’apprentissage, l’outil lui-même a ses limites :

• Résolution: Les simulations 3D sur une grille cartésienne sont coûteuses en calcul et nécessitent des grilles grossières, ce qui rend les résultats quantitatifs moins fiables.
• Turbulence: Le solveur est laminaire ; il lui manque un modèle de turbulence (comme k−ϵ) requis pour les écoulements complexes ou à grande vitesse.
• Diffusion: Les schémas de différenciation en amont sont stables mais numériquement diffusifs, potentiellement « étalés » sur les détails fins de l’écoulement.

Où aller à partir d’ici ? Ce projet sert de point de départ. Les améliorations futures pourraient inclure la mise en œuvre de schémas d’advection d’ordre supérieur (comme WENO), l’ajout d’une modélisation de la turbulence ou le passage à des méthodes de volumes finis (comme OpenFOAM) pour une meilleure gestion du maillage autour de géométries complexes. Il existe de nombreuses techniques astucieuses pour contourner la multitude de scénarios que vous souhaiterez peut-être mettre en œuvre. Ce n’est qu’un premier pas vers une véritable compréhension des CFD !