On s’intéresse à la résolution du système d’équations de Navier-Stokes
tridimensionnelles monophasiques, à une pression, instationnaires, en
incompressible ou faiblement dilatable, basées sur une discrétisation
temporelle de type Euler implicite d’ordre 1 ou Crank-Nicolson d’ordre 2 et sur
une discrétisation spatiale par volumes finis colocalisée.
Dans ce sous-programme sont calculées, à un pas de temps donné, les
variables vitesse et pression de ce problème en procédant en
deux étapes issues d’une décomposition des opérateurs (méthode à
pas fractionnaires).
Les variables sont donc supposées connues à
l’instant tn et on cherche à les déterminer à l’instant1 tn+1. Soit Δ tn =tn+1− tn le pas de temps associé. Dans un premier temps, on réalise l’étape de
prédiction de la vitesse en résolvant l’équation de quantité de
mouvement avec une pression explicite. Suit l’étape de correction de la
pression (ou projection de la vitesse) qui permet d’obtenir un champ de vitesse à divergence nulle.
Les équations en continu sont donc :
avec ρ la masse volumique, u le champ de vitesse,
[ TS−K u ] les autres termes sources (K est un
tenseur diagonal positif par définition), σ le tenseur
de contraintes, τ le tenseur des contraintes visqueuses, µ la
viscosité dynamique (moléculaire et éventuellement turbulente), κ
la viscosité de
volume (usuellement nulle et négligée dans le code et dans la suite du document,
sauf en compressible),
D le tenseur taux de déformation2, Γ le terme source de masse.
On rappelle la définition des notations employées3 :
et donc :
Dans le cas de la prise en compte d’une masse volumique variable, l’équation de continuité s’écrit :
| + div |
| = Γ |
Cette équation n’est pas prise en compte dans l’étape de projection (on continue à résoudre seulement div( ρ u) = Γ) alors que le terme ∂ ρ/∂ t apparaît lors de l’étape de prédiction de la vitesse dans le sous-programme preduv . Si ce terme joue un rôle sensible, l’algorithme compressible de Code Saturne(qui résout l’équation complète) est alors sans doute plus adapté.
On pose :
8cm[height=4cm]./graphics/facette8cm[height=4cm]./graphics/facebord
Figure II.0.1: Définition des différentes entités géométriques pour les faces internes (gauche) et de bord (droite).
Une des méthodes permettant de résoudre numériquement les équations de Navier-Stokes est de décomposer les opérateurs s’y rattachant en opérateurs moins complexes (qui peuvent être traités à l’aide d’algorithmes efficaces), moyennant des sous-pas intermédiaires dans un même pas de temps. Ici, deux sous-pas sont réalisés : le premier reprend les parties convective, diffusive et termes sources de l’équation de quantité de mouvement et constitue l’étape dite de prédiction de la vitesse, le second traite l’équation de continuité et est désigné comme l’étape de correction de pression ou de projection de la vitesse.
La discrétisation en temps se fait en appliquant à la variable résolue un θ-schéma au temps n+θ s’inspirant de la démarche utilisée pour l’équation de transport d’un scalaire4.
La vitesse au temps n+1 n’étant disponible qu’après l’étape de projection, c’est ici une vitesse prédite au temps n+1 que l’on utilise pour interpoler :
Avec5 :
Le champ de vitesse un+1 prédit est alors obtenu par
:
• Linéarisation partielle de l’opérateur de convection engendrant
un découplage des composantes de la vitesse.
• Explicitation de la pression.
• Explicitation ou extrapolation des grandeurs physiques (i.e
ρ, µ, Cp ...) et du flux de masse.
• Explicitation ou extrapolation des termes sources explicites au
temps n+θS tels que les contributions du gradient transposé, de la
viscosité secondaire, de la partie extradiagonale des pertes de charges, de la
masse injectée Γ ui, ...
• Les termes sources implicites linéaires par rapport à la vitesse
(termes sources utilisateurs implicite, partie diagonale des pertes de charge
Kd u, sources de masse Γ u, etc...)
sont supposés pris au temps n+θ et sont rattachés au tenseur
K6.
Par souci de clarté, on suppose qu’en l’absence d’indication, les propriétes
physiques Φ = ρ, µ ... et le flux de masse (ρu) pris
respectivement aux instants n+θΦ et n+θF, où θΦ et
θF dépendent des schémas en temps spécifiquement utilisés pour ces
grandeurs7.
On résout donc le système suivant, après réécriture des termes instationnaires à l’aide de la conservation de la masse :
avec :
Le champ de vitesse prédit est a priori à divergence non nulle. La seconde
étape corrige la pression en imposant la nullité de la contrainte
stationnaire pour la vitesse prise à l’instant tn+1.
On résout donc :
où l’incrément de pression δ Pn+θ vaut :
Les quantités ρ et µ sont constantes lors de ces deux étapes. Leur variation éventuelle est effectuée au début du pas de temps suivant, après réactualisation des scalaires (température, fraction massique,...).
On intègre classiquement sur les volumes de contrôle Ωi (ou cellules) les équations discrétisées en temps.
Si l’on ne tient pas compte des termes de convection et de diffusion issus du θ-schéma, les termes sources volumiques explicites de l’équation (??) s’écrivent pour le système portant sur la quantité (un+1−un) :
Pour intégrer ces termes sur une cellule Ωi, on multiplie leur valeur locale au centre par le volume de la cellule.
Après décomposition de u à l’aide de la relation (??), l’intégration spatiale des parties convectives θ div(un+1⊗ (ρ u)) et (1−θ) div( un⊗ (ρ u)) conduit à une somme de flux numériques F ij calculés aux faces des cellules purement internes et de flux numériques F bik calculés aux faces de bord du domaine Ω. Soient Vois(i) l’ensemble des centres des cellules voisines de Ωi et γb(i) l’ensemble des centres des faces de bord de Ωi, on a :
en posant :
La valeur de F ij dépend du type de schéma numérique
choisi. Il en existe trois dans Code Saturne :
| • un schéma décentré amont d’ordre 1 (upwind) : | |||||||
| F ij((ρ u), u)= F ij amont((ρ u),u) | |||||||
| où : | u f,ij=u I si (ρ u) ij.S ij 0 | ||||||
| u f,ij=u J si (ρ u) ij.S ij < 0 , | |||||||
| • un schéma centré d’ordre 2: | |||||||
| F ij((ρ u), u)= F ijcentré((ρu),u) | |||||||
| avec : | u f,ij= αij u I′+(1−αij) u J′ et uK′ = uK+(grad u)K. KK′ pour K = I ou J | ||||||
| • un schéma décentré amont d’ordre 2 SOLU (Second Order Linear Upwind) : | |||||||
| F ij((ρ u), u)= F ij SOLU((ρ u),u) | |||||||
| |||||||
La valeur de F bik est calculée avec :
avec u bik valeur au bord donnée directement par les
conditions aux limites.
En centré, on écrit en réalité (égalité conservant l’ordre un en espace) :
pour des raisons de stabilité purement numérique.
Un test de pente (qui peut introduire des non linéarités dans l’opérateur de convection) permet de basculer entre un schéma d’ordre deux et le schéma décentré amont d’ordre un. De plus, en mode standard, on utilise en tout point une valeur de u f,ij issue d’une moyenne barycentrique entre la valeur décentrée amont et la valeur d’ordre 2 (blending, spécifié par l’utilisateur).
De même, les parties diffusives θ div(µ grad un+1) et (1−θ) div(µ grad un) s’écrivent :
avec :
et :
en conservant les notations précédentes et avec
u bik la valeur au bord donnée directement par les conditions aux limites.
La viscosité µ ij à la face est calculée à l’aide des valeurs
aux centres selon une fonction f donnée :
qui est, soit une moyenne arithmétique :
soit une moyenne géométrique :
et la viscosité µ bik est égale à :
On introduit en outre, pour une utilisation ultérieure,
les notations suivantes :
qui correspondent chacune à une valeur non reconstruite aux faces internes et de bord.
Le système (??) pouvant comporter des non
linéarités dues au recours au test de pente ou pouvant conduire via la
reconstruction du gradient (cellule) à une matrice quasiment pleine en
présence de non orthogonalités, on le résout de manière itérative avec
la suite (un+1,k)k∈ IN définie par :
ce qui définit également la suite (δuk+1)
k∈ IN.
Les deux opérateurs E et EM ont pour expression respective :
avec :
De plus, on suppose que cette suite (un+1,k)k∈ IN
converge vers un+1.
Les conditions aux limites associées aux opérateurs E et EM du système (??) sont celles portant sur la vitesse u. Elles sont de type Dirichlet homogène ou de type Neumann homogène sur δu si u a une condition de type Dirichlet ou de type Neumann respectivement. Elles sont mixtes dans le cas d’une condition de symétrie sur une face en biais par rapport aux axes.
Les deux premières sommes de type (∑k∈ γb(i)),
i.e. comportant les termes en F bik((ρ u),u) et D bik(µ,u), utilisent les conditions aux limites de la vitesse.
Le terme volumique :
de terme de bord associé :
a un traitement particulier. En effet, pour une cellule Ωi jouxtant le
bord, la contribution du premier terme (en gradient transposé) est annulée, aucune condition à la limite correcte ne lui étant attribuée pour le moment.
L’opérateur EM approche E (aucun terme n’est reconstruit et la partie convective est traitée systématiquement en schéma décentré amont). Ceci peut générer des imprécisions numériques non négligeables si la suite (un+1,k)k∈ IN n’a pas convergé.
En prenant la divergence de la première équation du système
(??), on obtient :
En utilisant la contrainte stationnaire div(ρu)n+1=Γ, on a donc :
soit :
et :
En intégrant sur une cellule :
et :
On utilise le même formalisme que précédemment pour l’intégration du
terme de diffusion de l’étape de prédiction. Les conditions aux limites sont
de type Dirichlet homogène ou de type Neumann homogène sur δ P si P a une condition de
type Dirichlet ou de type Neumann respectivement.
La discrétisation de ∑j∈ Vois(i)[ (ρ
u)ijn+1. S ij] + ∑k∈ γb(i) [(ρ u) bikn+1 . S bik] est particulière. Le
choix suivant noté [ ]Init, pour une cellule ne touchant pas
le bord par exemple, est insatisfaisant avec la discrétisation et le schéma utilisés ici, en particulier avec l’équation
(??) :
Tout comme pour le calcul du flux numérique en centré, on
utilise de façon licite l’approximation suivante :
mais ce n’est pas elle qui pose problème.
En fait, [(ρu) ijn+1]Init contient le terme
G cel,ij n, hérité de l’étape de prédiction, qui vaut :
Or, sur un maillage orthogonal régulier cartésien, à partir d’une vitesse un+1=0 et d’une pression PI n−1+θ=(−1)I, on obtient G cel,ij n−1+θ=0 d’où δ Pn+θ=0 : l’irrégularité initiale de pression ne peut donc jamais être corrigée.
Pour remédier à cela, on modifie l’écriture [ ]Init
de (ρ
u) ijn+1 et de (ρ u) bikn+1 en
adoptant la valeur [ ]Corr :
et,
pour les conditions aux limites d’entrée, de symétrie (quelconque) ou de
paroi :
pour les conditions aux limites de sortie :
β est appelé coefficient d’Arakawa. Lorsqu’il vaut 1 (valeur par défaut), il s’agit du filtre Rhie & Chow.
On peut généraliser cette démarche à d’autres termes sources du même type, par exemple pour le modèle Rij−ε.
On construit une suite (δ Pn+1,k)k∈ IN pour résoudre
l’équation (??), qui peut conduire via la
reconstruction du gradient cellule à une matrice quasiment pleine en
présence de non orthogonalités, définie par :
ce qui définit également la suite (δ(δ P)n+θ,k+1)k∈
IN.
Les opérateurs F et FM ont pour expression :
et :
respectivement.
Crelax est un coefficient de relaxation donné et fixé à 1 en standard.
On suppose que la suite (δ Pn+θ,k)k∈ IN
converge
vers δ Pn+θ.
Au fur et à mesure des itérations, le flux de masse est mis à jour, en utilisant δ(δ P). À convergence, le flux de masse actualisé obtenu est :
et on calcule le nouveau champ de vitesse au centre des cellules grâce à l’égalité :
Un traitement spécifique permet d’assurer que la conservation de la masse portant sur le bilan des flux de masse aux faces est toujours parfaitement vérifiée à l’issue de l’étape de correction, que la suite (δ Pn+θ,k+1)k∈ IN ait ou non atteint la convergence. En effet, δ Pn+θ,kfin+1, dernier terme évalué de la suite, est donné par :
Au lieu de réactualiser classiquement le flux de masse, on le calcule comme
suit :
et :
ce qui conduit bien à :
soit :
et donc :
On se reportera aux sections relatives aux sous-programmes Preduv (prédiction des vitesses) et Resolp (correction de pression). Pour la reconstruction des vitesses par moindres carrés à partir des flux de masse aux faces, on pourra voir Recvmc .