Dans ce sous-programme, on effectue l’étape de prédiction de la vitesse
u. Ceci consiste à résoudre l’équation de quantité de
mouvement (Q.D.M.) en traitant la pression p de manière explicite. La solution en vitesse-pression est obtenue après une étape de correction sur la pression
effectuée dans le sous-programme resolp
, en utilisant la loi de conservation de la masse :
où Γ est le terme source de masse1.
L’équation de conservation de la quantité de mouvement moyenne obtenue par application
du théorème fondamental de la dynamique est :
où :
avec pour les écoulements newtoniens, la relation linéaire suivante :
σ représente le tenseur de contraintes, τ le tenseur
des contraintes visqueuses, µ la viscosité dynamique (moléculaire et
éventuellement turbulente), D
le tenseur taux de déformation2,
R le tenseur de Reynolds qui apparaît lors de l’application de
l’opérateur moyenne à l’équation instantanée, S les termes
sources.
λ est le second coefficient de viscosité. Il est relié à la viscosité de
volume κ par la relation
Quand l’hypothèse de Stokes est vérifiée, la viscosité de volume κ est
nulle, soit 3λ+2µ=0. Cette hypothèse est implicite dans le code et
dans le reste du doument, sauf en compressible.
L’équation de conservation de la quantité de mouvement s’écrit finalement
:
avec p définissant l’écart à la pression hydrostatique de référence (la
pression hydrostatique réelle étant calculée avec la masse volumique ρ
et non ρ 0) :
(X étant le vecteur de composantes x, y et z).
µt, K pdc, u i représentent respectivement
la viscosité dynamique turbulente, le tenseur des pertes de charge et la valeur de la variable associée à la source de
masse.
La divergence du tenseur des contraintes de Reynolds s’écrit :
Le terme source de masse fait intervenir la vitesse locale u et
aussi une vitesse u i associée à la masse injectée (ou retirée).
Lorsque Γ<0, on ôte de la masse au système et on a donc
u i = u. Le terme est nul (i.e. Γ
(u i−u)= 0 ). Quand Γ>0, on a un
terme non nul si u i ≠ u.
Dans ce sous-programme, tous les termes intervenant dans
l’équation de conservation de la
quantité de mouvement, excepté les termes de convection et diffusion, sont
calculés et transmis au sous-programme codits
qui résout l’équation complète
(convection-diffusion avec termes sources).
Le terme convectif en div(u ⊗ ρ u)
introduit une non linéarité et un couplage des composantes de la vitesse
u dans l’équation (??). Une linéarisation et un découplage
des trois composantes de la
vitesse sont réalisés lors de la discrétisation de cette étape de
prédiction.
En effet, soit :
La contribution exacte du terme convectif à prendre en compte dans cette
étape de prédiction serait :
Les deux derniers termes de l’expression (??) sont a priori négligés
de manière à obtenir un système en vitesse qui soit découplé et donc,
éviter l’inversion d’une matrice pouvant être de très grande taille. Ces
deux termes peuvent néanmoins être pris en compte de manière plus ou moins
approchée par extrapolation explicite du flux de masse en n+θF (pour le
terme couplant linéaire provenant de la convection de un par δ
u) et utilisation d’un point-fixe par sous itération sur le sous
programme navsto
(pour le terme non-linéaire, en spécifiant NTERUP>1).
L’équation (??) est discrétisée au temps n+θ à l’aide d’un
θ-schéma, et le tenseur des pertes de charges décomposé en une partie
diagonale Kd et une extradiagonale Ke (soit
Kpdc=Kd+Ke).
• La pression est supposée connue en n−1+θ (décalage temporel
pression-vitesse) et le gradient naturellement calculé à cet instant.
• Les termes sources de viscosité secondaire, de gradient transposé,
ceux provenant du modèle de turbulence3, ρ K e u, (ρ −ρ0)
g ainsi que Ts exp et
Γ u i sont pris de manière explicite au temps n, ou
extrapolés suivant le schéma en temps choisi pour les propriétés physique et les
termes sources.
• Les termes sources u div (ρ u),
Γ u, Ts imp u et
−ρ K d u sont implicités est calculés à
l’instant n+θ.
• Le terme de diffusion div (µ tot grad u) est
implicité : la vitesse est prise à l’instant n+θ et la viscosité
explicitée ou extrapolée.
• Enfin, le terme de convection en div( u ⊗
(ρu) ) est implicité : la composante résolue de la vitesse est
prise en n+θ, et le flux de masse, explicité, ou extrapolé en
n+θF.
Par souci de clarté, on suppose, en l’absence d’indication, que les propriétes physiques Φ (ρ, µtot, ...) et le flux de masse (ρu) sont pris respectivement aux instants n+θΦ et n+θF, où θΦ et θF dépendent des schémas en temps spécifiquement utilisés pour ces grandeurs4.
La discrétisation temporelle de l’équation (??) s’écrit alors comme suit :
où, par souci de simplification, on a posé :
et :
Par analogie avec l’écriture du θ-schéma pour une variable scalaire,
un+θ est interpolée à partir de la vitesse
prédite un+1 de la manière suivante5 :
Avec :
L’équation (??) est alors réécrite sous la forme :
d’où l’équation résolue par le sous-programme codits
:
La méthode de discrétisation spatiale est développée dans le sous-programme codits
.
Dans le cas standard sans extrapolation, le terme
− Ts imp n’est ajouté à fs imp que s’il est positif afin de ne
pas affaiblir la dominance de la diagonale de la matrice à inverser.
Si une extrapolation est utilisée, par contre,
Ts imp est ajouté à fs imp quel que soit son signe. En effet, l’idée intuitive qui
consiste à prendre :
aboutit à une incohérence dans le traitement si Tsimp change de signe
entre deux pas de temps.
la partie diagonale K d du terme
de perte de charge est utilisée dans fs imp. En fait, pour être rigoureux,
il faudrait ne retenir que les contributions positives (point signalé dans le
sous-programme utilisateur associé uskpdc
). Cette prise en compte sera à améliorer.
Le terme θ Γn−θ div
(ρ u) ne pose pas de problème pour la
dominance de la diagonale de la matrice car il est exactement compensé par le
terme de convection (cf. covofi
).
L’équation de conservation de la quantité de mouvement est donc résolue de façon
découplée. Ainsi, l’intégration des différents termes a été effectuée afin de
traiter séparément l’équation obtenue pour chaque composante de la vitesse.
Dans le sous-programme preduv
, on calcule pour chaque
composante le second membre fsexp du système (??), les termes implicites du système (à l’exception des termes de convection-diffusion), et le terme de viscosité totale aux faces internes6 et de
bord. Ces termes sont alors transmis au sous-programme codits
qui construit et
résout le système complet obtenu pour chaque composante de la vitesse avec les termes de convection-diffusion.
Le résidu de normalisation pour la résolution du système en pression (resolp
) est calculé dans
preduv
. Il est défini par la norme de la
grandeur
| div(ρ |
| n+1+Δ t g |
| Pn−1+θ)−Γ |
intégrée sur chaque cellule IEL du maillage (Ωiel) soit, symboliquement, par la racine carrée de la somme sur les cellules du maillage de la quantité
Il représente le second membre du système qui porterait sur la pression si le gradient de pression n’était pas pris en compte lors de l’étape de prédiction des vitesses. On note que si l’on utilisait directement le second membre de l’équation portant sur l’incrément de pression, on obtiendrait, pour un calcul stationnaire mené à convergence, un résidu de normalisation tendant vers zéro, ce qui serait pénalisant et peu utile.
Au début de preduv , on ne dispose pas encore de un+1 et il n’est donc pas possible de calculer le résidu de normalisation en totalité. Cependant, le calcul du résidu complet à la fin de preduv n’est pas souhaitable non plus, car on devrait alors monopoliser un tableau de travail pour conserver le gradient de pression tout au long de preduv . Le calcul du résidu de normalisation est donc réalisé en deux fois.
La quantité ∫Ωieldiv(Δ t grad Pn−1+θ) −Γ dΩ est calculée au début de preduv et on y ajoute le complément ∫Ωieldiv(ρ un+1)dΩ à la fin de preduv .
On calcule donc tout d’abord le gradient de pression aux cellules à l’instant n−1+θ par un appel à grdcel . On utilise alors inimas pour évaluer Δ t S grad Pn−1+θ·n aux faces (de surface S et de normale n). Pour cela, en entrée de inimas , le tableau de travail TRAV contient Δ t/ρ grad Pn−1+θ ; en sortie, les tableaux VISCF et VISCB contiennent la valeur de Δ t S grad Pn−1+θ·n aux faces internes et de bord respectivement.
On utilise ensuite divmas qui place alors dans XNORMP la valeur de ∫Ωieldiv(Δ t grad Pn−1+θ) dΩ aux cellules à partir des tableaux VISCF et VISCB.
On ajoute enfin à XNORMP la contribution ∫Ωiel −Γn dΩ du terme source de masse.
On applique pour ρ u + Δ t grad P les conditions aux limites de la vitesse. Les conditions aux limites utilisées pour le gradient de pression (ou plutôt pour Δ t/ρ grad Pn−1+θ) pour le calcul de ∫Ωieldiv(Δ t grad Pn−1+θ) dΩ sont donc les conditions aux limites de la vitesse homogénéisées : ainsi, on suppose que dans la direction normale aux entrées et aux parois, le gradient de pression (ou plutôt Δ t/ρ grad Pn−1+θ) est nul et que dans la direction normale aux symétries et aux sorties, il reste inchangé.
De plus, pour gagner du temps calcul lors du passage par inimas , on se contente, sur les maillages non orthogonaux, d’une évaluation des valeurs aux faces à l’ordre 1 en espace (pas de reconstruction : NSWRP=1). En effet, on cherche à évaluer un simple résidu de normalisation global : la précision locale n’a donc pas d’intérêt.
Le calcul du résidu sera complété à la fin de preduv .
• Calcul en partie du résidu de normalisation pour l’étape de pression
Dans cette première étape on calcule dans le tableau XNORMP(NCELET) la
grandeur
| div(Δ t g |
| Pn−1+θ)−Γ |
intégrée sur chaque cellule IEL du maillage (Ωiel) soit, symboliquement,
| XNORMP(IEL)= |
|
|
On réalise cette opération en utilisant successivement inimas
(calcul aux faces dans VISCF et VISCB de Δ t grad Pn−1+θ à
partir du tableau de travail TRAV=Δ t/ρ
grad Pn−1+θ, assorti des conditions aux limites de vitesse
homogènes et sans reconstruction) et divmas
(calcul dans XNORMP
de l’intégrale sur les cellules). Par une simple boucle, on
ajoute ensuite la contribution du terme source de masse Γ.
Ce calcul est complété à la fin de preduv
.
• Calcul en partie du terme fs exp
Pour représenter le second membre correspondant à chaque composante de la
vitesse, on utilise les tableaux TRAV(IEL,DIR),
TRAVA(IEL,DIR) et PROPCE,
où IEL est le numéro de la cellule et DIR la direction (x, y,
z). Quatre cas sont à considérer suivant que les termes sources sont extrapolés
en n+θS, ou que l’on itère par un point fixe sur le système en
vitesse-pression (NTERUP>1).
• Si on extrapole les termes sources et que l’on itère sur navsto
• Sans itération sur navsto
, TRAVA est inutile et son contenu est directement stocké dans TRAV.
• Sans extrapolation des termes sources, PROPCE est inutile et son contenu est directement stocké dans TRAVA (ou dans TRAV si TRAVA est inutile).
Ainsi, sans extrapolation des termes sources, et sans itération sur navsto
, tout les termes sources vont directement dans TRAV.
• Calcul du terme de viscosité aux faces
(µ tot)ijSURFN/DIST
Le calcul du terme de viscosité totale aux faces est effectué par le
sous-programme viscfa
et stocké dans les tableaux VISCF et
VISCB pour les faces internes et faces de bord respectivement.
Lors de l’intégration des termes de convection-diffusion dans le sous-programme
codits
, on distingue les termes non reconstruits, intégrés dans la
matrice EM , de l’ensemble des termes (non reconstruits +
gradients de reconstruction) associés à l’opérateur E (non
linéaire)9.
De la même manière, on distingue la viscosité totale aux faces utilisée dans
E, tableaux VISCF et VISCB, de la viscosité totale aux
faces utilisée dans EM, tableaux VISCFI et
VISCBI.
Pour les modèles à viscosité turbulente et en LES, ces deux tableaux sont identiques et contiennent µt+µ.
Pour les modèles au second ordre, ils contiennent normalement µ, mais pour des simples raisons de stabilité
numérique, on peut choisir de mettre µt+µ dans la matrice (i.e. dans EM) en
conservant µ au second menbre (i.e. dans E). De par la résolution par incréments, cette
manipulation ne change pas le résultat. Cette option est activée par l’indicateur IRIJNU = 1
Si la vitesse n’est pas diffusée ( IDIFF(IUIPH) < 1), alors les termes VISCF et VISCB sont mis à zéro.
• Calcul du second membre complet, de fs imp et résolution de l’équation
Les équations d’évolution des composantes de la quantité de mouvement sont
résolues de façon découplée. On utilise, par conséquent, un seul tableau
ROVSDT pour représenter la partie diagonale de la matrice obtenue pour chaque composante de la vitesse.
Pour chaque composante de la vitesse :
| SMBR=(1−θS) PROPCE+TRAVA+TRAV |
| SMBR=TRAVA+TRAV |
avec RHOn le tenseur diagonal d’élément ρIELn,
Ω le tenseur diagonal d’élément |ΩIEL|, 1 le
vecteur de composantes toutes égales à 1.
Fin de la boucle sur les composantes de la vitesse.
• Fin du calcul du résidu de normalisation pour l’étape de pression
Comme indiqué précédemment, on peut maintenant compléter le calcul du
résidu de normalisation pour l’étape de pression de resolp
.
Le tableau XNORMP contient déjà ∫Ωieldiv(Δ t grad Pn−1+θ) −Γ dΩ . On lui ajoute donc ∫Ωieldiv(ρ un+1)dΩ.
Pour cela, on procède comme précédemment pour le calcul de ∫Ωieldiv(Δ t grad Pn−1+θ) dΩ. Un appel à inimas permet d’obtenir ρ S un+1·n aux faces à partir de un+1 connu aux cellules (tableau RTP). Les conditions aux limites pour inimas sont naturellement celles de la vitesse. Comme précédemment, on se contente pour gagner du temps calcul lors du passage par inimas , d’une évaluation des valeurs aux faces à l’ordre 1 en espace sur les maillages non orthogonaux (pas de reconstruction : NSWRP=1). On utilise ensuite divmas pour calculer aux cellules la divergence ∫Ωieldiv(ρ un+1)dΩ et l’ajouter directement à XNORMP.
Pour finir, le résidu de normalisation est déterminé et stocké dans RNORMP(IIPHAS) par un appel à prodsc
(qui réalise le calcul de la somme sur les cellules du
carré des valeurs de XNORMP et en prend la racine carrée).
On résume dans les tableaux (??), (??), (??)
et (??) les différentes contributions (hors convection-diffusion)
affectées à chacun des vecteurs TRAV, TRAVA, PROPCE et
ROVSDT à l’itération n. On différencie pour chacun des schémas en temps
appliqués aux les termes sources, deux cas suivant qu’un point fixe sur le
système en vitesse-pression est utilisé ou non (itération sur navsto
pour
NTERUP>1). En l’absence d’indication les propriétés physiques Φ
(ρ,µ,etc...) sont supposées prises au temps n+θΦ, et le flux
de masse ( ρ u) pris au temps n+θF, où θΦ
et θF dépendent des schémas en temps spécifiquement utilisés pour ces
grandeurs (cf. introd
).
Les termes figurant dans ces tableaux sont écrits tels qu’ils ont été implantés
dans le code, d’où l’origine de certaines différences par rapport à l’écriture
adoptée dans l’équation (??).
Par souci de simplification, on pose :
• NTERUP = 1 : SMBRn=(1−θS) PROPCEn+TRAVn
• NTERUP > 1 (sous-itération k) : SMBRn=(1−θS) PROPCEn+TRAVAn+TRAVn
• NTERUP = 1 : SMBRn=TRAVn
• NTERUP > 1 (sous-itération k) : SMBRn=TRAVAn+TRAVn
• Prise en compte du terme grad (ρ k) pour les modèles à viscosité turbulente
Pour les modèles à viscosité turbulente, on calcule ρ grad k au lieu de grad (ρ k). Cette
approximation, historique, provient du fait que les conditions aux limites de
ρ k ne sont pas directement accessibles, contrairement à celles de
k.
• Prise en compte de la diagonale de K pdc
Actuellement, dans le sous-programme utilisateur uskpdc
, une mise en
garde explicite est écrite, mais en commentaire.
La partie diagonale K d du tenseur de pertes de charge
K pdc peut donc intervenir systématiquement dans le calcul du
coefficient fs imp, que sa contribution K d soit positive ou non, si
l’utilisateur n’y prend garde. Un test de positivité sur les éléments de
K d assurant une prise en compte correcte (contribution renforçant réellement la diagonale de la matrice globale) devrait être implanté.
• Écriture de EM
Dans la résolution procédant par incréments, il n’est pas indispensable
à convergence
que la viscosité utilisée pour l’écriture de l’opérateur E soit la même que celle prise en
compte dans EM, matrice du système en incréments. Ainsi, en Rij−ε, la
viscosité totale utilisée dans EM contient la viscosité moléculaire mais aussi la
viscosité turbulente si l’on choisit l’option IRIJNU = 1 , alors que
dans E intervient seule la viscosité moléculaire. Cet ajout de
la viscosité turbulente qui n’a pas de raison d’apparaître en
Rij−ε, a été hérité des pratiques mises en œuvre dans
ESTET et N3S-EF pour renforcer la stabilité (lissage éventuel de l’incrément). Mais, ce
n’est peut être pas le seul effet produit. En outre, cette pratique n’a pas aujourd’hui montré son absolue
nécessité dans Code Saturne. Par conséquent, une étude approfondie serait intéressante.
• Résidu de normalisation de l’étape de pression
On pourra vérifier le calcul du résidu de normalisation et en
particulier l’utilisation des conditions aux limites de vitesse.
• Calcul des pertes de charges
Avec extrapolation des termes sources on a :
Il serait aussi envisageable d’utiliser :
