Équations de Navier-Stokes et simulation numérique

Algorithme SIMPLE — Introduction à FLUENT

Author

École des Mines

1 Rappel et mise en contexte

1.1 Le cas d’étude du TD

Dans les séances 3 à 6, on étudie un injecteur air–éthane :

  • géométrie simplifiée en 2D cartésienne
  • air injecté en amont
  • éthane injecté à 45° à l’extrados du coude
  • régime permanent, laminaire, isotherme
  • sans turbulence et sans réaction chimique

Objectif physique : analyser le mélange des deux gaz par diffusion et convection.

1.2 Pourquoi ce cours avant la séance 4 ?

Contexte

La séance 4 du TD consiste à mettre en données puis à calculer l’écoulement dans FLUENT.

Jusqu’ici, FLUENT a été utilisé comme un outil de calcul.
On cherche maintenant à comprendre le problème numérique qu’il résout.

WarningDifficulté centrale

Pour ce cas, il faut déterminer simultanément :

  • le champ de vitesse \(\vv{u}\)
  • le champ de pression \(p\)
  • avec la contrainte de continuité \[ \nabla\cdot\vv{u}=0 \]

Or la pression n’est pas obtenue par une équation de transport indépendante.

ImportantQuestion du cours

Comment obtenir numériquement un champ de vitesse et un champ de pression compatibles ?

Réponse proposée :
présenter l’algorithme SIMPLE, méthode historique de couplage pression–vitesse, à la base des solveurs de type pressure-based.

1.3 La référence fondatrice

NoteSuhas V. Patankar
  • Doctorant de Brian Spalding à Imperial College (PhD 1967)
  • Développe la méthode des volumes finis appliquée aux écoulements
  • Propose avec Spalding l’algorithme SIMPLE (1972)

Ouvrage de référence

Numerical Heat Transfer and Fluid Flow
Patankar, 1980 — ISBN 0-89116-522-3

Cet ouvrage reste une référence majeure en mécanique des fluides numérique.

2 Les équations gouvernantes

2.1 L’équation de transport générale

NoteIdée clé : toutes les variables obéissent au même schéma

Soit \(\phi\) une quantité transportée (vitesse, température, concentration…)

\[ \boxed{ \underbrace{\pd{(\rho\phi)}{t}}_{\text{accumulation}} + \underbrace{\nabla\cdot(\rho\,\vv{u}\,\phi)}_{\text{convection}} = \underbrace{\nabla\cdot(\Gamma_\phi\,\nabla\phi)}_{\text{diffusion}} + \underbrace{S_\phi}_{\text{source}} } \]

Conservation de masse : \(\phi=1\), \(\Gamma_\phi=0\), \(S_\phi=0\) \(\;\Rightarrow\;\) \(\pd{\rho}{t} + \nabla\cdot(\rho\vv{u}) = 0\)

2.2 Les équations de Navier-Stokes

NoteQuantité de mouvement : \(\phi = u_i\)

\[ \pd{(\rho u_i)}{t} + \pd{(\rho u_j u_i)}{x_j} = -\pd{p}{x_i} + \pd{\tau_{ij}}{x_j} + \rho g_i \]

Pour un fluide newtonien incompressible : \[ \tau_{ij}=\mu\!\left(\pd{u_i}{x_j}+\pd{u_j}{x_i}\right) \qquad\Rightarrow\qquad \pd{\tau_{ij}}{x_j}=\mu\nabla^2 u_i \]

Inconnues du problème (3D)

Variable Rôle
\(u,\,v,\,w\) composantes de vitesse
\(p\) pression
Total 4 inconnues

3 équations de quantité de mouvement + 1 équation de continuité.

WarningDifficultés principales
  1. Non-linéarité des termes convectifs
  2. Couplage entre vitesse et pression
  3. Absence d’équation directe pour la pression
  4. Nécessité d’une résolution itérative

2.3 Couplage pression–vitesse

NoteÉquations de quantité de mouvement

\[ \rho\pd{\vv{u}}{t} + \rho(\vv{u}\cdot\nabla)\vv{u} = -\nabla p + \mu\nabla^2\vv{u} \]

  • inconnues : \(\vv{u}\) et \(p\)
  • la pression intervient par le terme \(\nabla p\)
NoteÉquation de continuité

\[ \nabla\cdot\vv{u} = 0 \]

  • contrainte d’incompressibilité
  • la pression n’y apparaît pas explicitement
WarningConséquence
  • Il n’existe pas d’équation explicite donnant directement le champ de pression.
  • Le champ de pression doit être tel que le champ de vitesse issu des équations de quantité de mouvement vérifie la continuité.

\(\Rightarrow\) le problème numérique consiste à assurer la cohérence entre pression et vitesse.

2.4 Couplage pression–vitesse : idée essentielle

Problème fondamental

Si le champ de pression est connu, les équations de quantité de mouvement permettent de calculer la vitesse.
En revanche, la détermination du champ de pression compatible avec la continuité n’est pas directe.

Interprétation

La pression joue le rôle de variable de contrainte assurant la satisfaction de la condition d’incompressibilité.

Conséquence numérique

La résolution nécessite une méthode spécifique de couplage pression–vitesse, telle que SIMPLE.

3 Du maillage au solveur

3.1 Discrétisation : la méthode des volumes finis

NotePrincipe : intégration sur un volume de contrôle \(V_P\)

\[\int_{V_P}\pd{(\rho\phi)}{t}\,dV + \oint_S \rho\vv{u}\phi\cdot\vv{n}\,dS = \ldots\]

  • formulation conservative par construction
  • adaptée aux géométries complexes
  • utilisée dans les solveurs industriels (FLUENT, OpenFOAM)

3.2 Maillage colocalisé : instabilité en damier

WarningOrigine

Si \(p\) et \(\vv{u}\) sont stockés au même nœud :

\[ \frac{\partial p}{\partial x}\bigg|_i = \frac{p_{i+1} - p_{i-1}}{2\Delta x} \]

Le point \(i\) n’intervient pas dans son propre gradient.

NoteConséquence
  • découplage entre nœuds pairs et impairs
  • existence de champs de pression oscillants non détectés

\(\Rightarrow\) solution numérique non physique

3.3 Problème du maillage colocalisé : l’instabilité en damier

WarningStockage des variables

Si \(p\) et \(\vv{u}\) sont stockés au même nœud (maillage colocalisé), le gradient de pression centré en \(i\) s’écrit : \[\frac{\partial p}{\partial x}\bigg|_i = \frac{p_{i+1} - p_{i-1}}{2\Delta x}\]

Le noeud \(i\) n’intervient pas dans l’expression de son propre gradient.

WarningConséquence

Les nœuds pairs et impairs se découplent et forment deux sous-grilles indépendantes.

Un champ de pression oscillant (haut–bas–haut…) constitue alors une solution numérique admissible mais non physique, que le solveur ne détecte pas.

→ perte d’unicité et dégradation de la solution.

3.4 La solution : le maillage décalé (staggered grid, Harlow & Welch 1965)

TipPrincipe
  • pression stockée au centre des cellules
  • \(u\) sur les faces verticales
  • \(v\) sur les faces horizontales

Le gradient pour \(u_{i+\frac{1}{2}}\) devient : \[\frac{\partial p}{\partial x}\bigg|_{i+\frac{1}{2}} = \frac{p_{i+1} - p_i}{\Delta x}\]

Le calcul fait intervenir des nœuds adjacents, ce qui rétablit le couplage.

NoteInterprétation

Toute oscillation du champ de pression génère un gradient non nul et est donc détectée par le schéma numérique.

3.5 La solution : le maillage décalé (staggered grid, Harlow & Welch 1965)

NoteComparaison
Maillage Gradient Instabilité damier
Colocalisé \((p_{i+1}-p_{i-1})/2\Delta x\) non détectée
Décalé \((p_{i+1}-p_i)/\Delta x\) détectée et corrigée

Le calcul du gradient de pression constitue le mécanisme central du couplage entre pression et vitesse.
Le choix du maillage est donc déterminant pour la stabilité et la qualité de la solution.

3.6 Équation de quantité de mouvement discrétisée

NoteVolume de contrôle centré sur la face \(e\) (pour \(u_e\))

\[a_e\, u_e = \sum_{nb} a_{nb}\, u_{nb} + (p_P - p_E)A_e + b_e\]

avec \(a_e\), \(a_{nb}\) : contributions convectives et diffusive

NoteDirection \(x\) : \(u_e\)

\[a_e\, u^*_e = \sum_{nb} a_{nb}\, u^*_{nb} + (p^*_P - p^*_E)A_e + b_e\]

NoteDirection \(y\) : \(v_n\)

\[a_n\, v^*_n = \sum_{nb} a_{nb}\, v^*_{nb} + (p^*_P - p^*_N)A_n + b_n\]

WarningPoint clé

Si la pression n’est pas connue, les composantes de vitesse ne peuvent pas être déterminées.

Une équation supplémentaire est nécessaire pour fermer le système : elle est obtenue à partir de la continuité \(\nabla\cdot\vv{u}=0\).

4 L’algorithme SIMPLE

4.1 Principe : estimation et correction itérative

SIMPLE = Semi-Implicit Method for Pressure-Linked Equations (Patankar & Spalding, 1972)

NotePrincipe

On considère un champ de pression estimé \(p^*\), à partir duquel on calcule un champ de vitesse \(u^*\), \(v^*\) qui ne satisfait pas, en général, la continuité.

Le champ de pression et les vitesses sont ensuite corrigés de manière itérative jusqu’à satisfaction de la conservation de la masse.

TipInterprétation

La condition de continuité joue le rôle de critère de cohérence du champ de vitesse calculé.

L’algorithme ajuste le champ de pression de façon à rendre ce champ de vitesse compatible avec l’incompressibilité.

4.2 Les équations de correction

NoteÉtape 1 — Résolution avec \(p^*\)

On résout les équations de quantité de mouvement avec la pression estimée : \(a_e\, u^*_e = \sum_{nb} a_{nb}\, u^*_{nb} + (p^*_P - p^*_E)A_e + b_e\)

NoteÉtape 2 — Décomposition et linéarisation

\[ p = p^* + p', \qquad u = u^* + u' \]

On écrit l’équation de quantité de mouvement pour la correction de vitesse \(u'\).

Cette équation contient en principe :

  • un terme local \(a_e u'_e\),
  • des contributions des mailles voisines \(\sum a_{nb} u'_{nb}\),
  • un terme de pression.

Dans l’approximation SIMPLE, on néglige les contributions voisines \(\sum a_{nb} u'_{nb}\).

→ On conserve uniquement une relation locale entre \(u'\) et \(p'\).

On obtient alors :

\[ u'_e = \frac{A_e}{a_e}(p'_P - p'_E) \qquad\Rightarrow\qquad u_e = u^*_e + \frac{A_e}{a_e}(p'_P - p'_E) \]

4.3 Les équations de correction

NoteÉtape 3 — Équation de correction de pression

La vitesse est exprimée comme : \[ u = u^* + u' \]

On injecte cette expression dans la continuité : \[ \nabla\cdot\vv{u} = 0 \]

  • le champ \(u^*\) est connu → sa divergence apparaît comme un terme source explicite ;
  • la correction \(u'\) dépend de \(p'\) → elle fournit les inconnues.

On obtient ainsi une équation de type Poisson pour \(p'\) :

\[ a_P\,p'_P = a_E\,p'_E + a_W\,p'_W + a_N\,p'_N + a_S\,p'_S + b_P \]

avec : \[ b_P \propto \nabla\cdot\vv{u}^* \]

\(b_P\) mesure le défaut local de continuité du champ provisoire.

WarningSous-relaxation

La sous-relaxation est nécessaire pour stabiliser l’algorithme. L’approximation introduite dans l’évaluation de \(u'\) tend sinon à produire une correction excessive.

\(p^{\text{new}} = p^* + \alpha_p\,p' \qquad \alpha_p \approx 0{,}3\)

4.4 Variantes et paramètres dans FLUENT

TipVariantes de SIMPLE
Méthode Caractéristique
SIMPLE Formulation de référence, robuste
SIMPLEC Prise en compte plus complète des termes voisins
SIMPLER Amélioration de l’estimation initiale de \(p^*\)
PISO Corrections successives, particulièrement adaptées au transitoire
WarningParamètres usuels dans FLUENT
Paramètre Valeur typique Rôle
\(\alpha_p\) 0,3 Sous-relaxation sur la pression
\(\alpha_u\) 0,7 Sous-relaxation sur la vitesse
Résidu \(10^{-4}\)\(10^{-6}\) Critère de convergence

Si la convergence est lente, on peut diminuer \(\alpha_p\).

Si les résidus oscillent, il est également possible de diminuer \(\alpha_u\).

5 Les solveurs de FLUENT

5.1 Deux familles de solveurs

Solveur à base de pression
(Pressure-based)

  • Algorithmes : SIMPLE, SIMPLEC, PISO
  • Traitement du couplage pression–vitesse
  • Adapté en priorité aux écoulements incompressibles ou faiblement compressibles

Exemples : ventilation, pompe, aile subsonique

Solveur à base de densité
(Density-based)

  • Traitement plus direct des écoulements compressibles
  • Couplage fort avec l’équation d’état
  • Adapté en priorité aux écoulements où les variations de densité sont importantes

Exemples : tuyère, onde de choc, écoulement supersonique

Remarque : le nombre de Mach donne un ordre de grandeur utile, mais ne constitue pas un critère absolu de choix.

6 Transports couplés

6.1 Les équations sont-elles indépendantes ?

On considère :

  • Navier–Stokes → \(\vv{u}, p\)
  • Énergie → \(T\)
  • Espèces → \(Y_k\)

6.1.1 Question

Peut-on résoudre ces équations séparément ?

6.1.2 Observation

Les propriétés physiques interviennent dans les équations :

  • viscosité \(\mu\), masse volumique \(\rho\), diffusivité \(D\)

Or :

  • \(\mu = f(T, Y_k)\) , \(\rho = f(T, Y_k)\)

👉 Les équations sont liées entre elles

6.2 Notion de transports couplés

Une variation de température ou de composition :

  • modifie les propriétés du fluide
  • modifie l’écoulement

6.2.1 Exemples

  • \(T \uparrow\)\(\mu\) varie → modification des contraintes visqueuses
  • \(T \uparrow\)\(\rho\) varie → apparition de flottabilité
  • \(Y_k\) varie → \(\rho\) varie → modification de la dynamique

6.2.2 Conclusion

Les équations :

  • ne sont pas indépendantes
  • interagissent via les propriétés

👉 On parle de transports couplés

7 Vers le TD injecteur

7.1 Ce que vous allez utiliser dans FLUENT

Pour la séance 4, le cadre de calcul est le suivant :

  • solveur pressure-based
  • algorithme de couplage SIMPLE
  • régime stationnaire
  • modèle laminaire
  • module Species Transport
  • mélange air–éthane
  • schéma convectif upwind du premier ordre

7.1.1 Sens physique de ces choix

  • laminaire : le régime d’écoulement doit être vérifié par le nombre de Reynolds
  • isotherme : pas de couplage avec l’énergie à ce stade
  • sans combustion : on isole le phénomène de mélange
  • SIMPLE : on traite d’abord correctement le couplage entre vitesse et pression

7.2 Ce que le TD devra montrer

À l’issue de la séance 4, vous devrez être capables de :

  • vérifier que l’écoulement est bien laminaire
  • obtenir un calcul convergé
  • visualiser le champ de vitesse
  • visualiser la fraction massique d’éthane
  • vérifier les bilans de matière
  • tracer des profils de composition dans le tronçon de sortie
  • analyser l’effet du coefficient de diffusion

7.2.1 Question physique finale

Comment le mélange air–éthane évolue-t-il le long du coude, et comment l’épaisseur de diffusion se développe-t-elle vers la sortie ?