Le but de ce sous-programme est principalement de calculer le flux de masse aux faces. Il prend une variable vectorielle associée au centre des cellules (généralement la vitesse), la projette aux faces en la multipliant par la masse volumique, et la multiplie scalairement par le vecteur surface. Plus généralement, inimas est aussi appelé comme première étape du calcul d’une divergence (terme en div(ρR) en Rij−ε, filtre Rhie & Chow, ...).
La figure II.0.1 rappelle les diverses définitions géométriques pour les faces internes et les faces de bord. On notera α=FJ′/I′J′ (défini aux faces internes uniquement).
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).
On ne connaît pas la masse volumique à la face, cette dernière doit donc aussi être interpolée. On utilise la discrétisation suivante :
La partie en α(ρI uI)
+(1−α)(ρJ uJ) correspondant en fait à
(ρu)O. Le gradient en O est calculé par interpolation :
grad (ρu)O=
1/2[grad (ρu)I+grad (ρu)J]. La
valeur 1/2 s’est imposée de manière heuristique au
fil des tests
comme apportant plus de stabilité à l’algorithme global qu’une interpolation
faisant intervenir α. L’erreur commise sur ρu est en
O(h2).
Le traitement des faces de bord est nécessaire pour y calculer le flux de masse, bien sûr, mais aussi pour obtenir des conditions aux limites pour le calcul du grad (ρ u) utilisé pour les faces internes.
Pour les faces de bord, on connaît la valeur de ρF, qui est stockée
dans la variable ROMB. De plus, les conditions aux limites pour u
sont données par des coefficients A et B tels que :
(k∈{1,2,3} est la composante de la vitesse, l’erreur est en O(Bkh))
On a donc à l’ordre 1 :
Mais pour utiliser cette formule, il faudrait calculer grad (u) (trois
appels à GRDCEL
), alors qu’on a déjà calculé
grad (ρu) pour les faces internes. Le surcoût en temps serait alors
important. On réécrit donc :
Pour calculer les gradients de ρu, il faudrait donc en théorie
utiliser les coefficients de conditions aux limites équivalents :
Ãk = ρF Ak
Bk = BkρF/ρI′
Ceci paraît délicat, à cause du terme en
ρF/ρI′, et en particulier à l’erreur
que l’on peut commettre sur ρI′ si la reconstruction des gradients
est imparfaite (sur des maillages fortement non orthogonaux par exemple).
On réécrit donc l’équation
(??) sous la forme suivante :
Pour le calcul du flux de masse au bord, on va faire deux approximations. Pour
le deuxième terme, on va supposer ρI′≈ρI (ce qui
conduit à une erreur en O(Bkh) sur ρu si
ρI′≠ ρI). Pour le
troisième terme, on va supposer ρI′≈ρF. Cette
dernière approximation est plus forte, mais elle n’intervient que dans la
reconstruction des non-orthogonalités ; l’erreur finale reste donc faible
(erreur en O(Bkh2) sur ρu si
ρI′≠ ρF).
Et au final, le flux de masse au bord est calculé par :
Pour le calcul des gradients, on repart de l’équation (??), sur
laquelle on fait l’hypothèse que ρI′≈ρF. Encore une
fois, cette hypothèse peut être assez forte, mais les gradients obtenus ne
sont utilisés que pour des reconstructions de non-orthogonalités ; l’erreur
finale reste donc là encore assez faible.
Au final, les gradients sont calculés à partir de la formule suivante :
ce qui revient à utiliser les conditions aux limites suivantes pour
ρ u:
Ãk = ρF Ak
Bk = Bk
Dans la plupart des cas, les approximations effectuées n’engendrent aucune
erreur. En effet :
- dans le cas d’une entrée on a généralement Bk=0, avec un flux de
masse imposé par la condition à la limite.
- dans le cas d’une sortie, on a généralement flux nul sur les scalaires
donc sur ρ, soit ρF=ρI′=ρI.
- dans le cas d’une paroi, on a généralement Bk=0 et le flux de masse
est imposé nul.
- dans le cas d’une symétrie, on a généralement
ρF=ρI′=ρI et le flux de masse est imposé nul.
Pour sentir un effet de ces approximations, il faudrait par exemple une paroi
glissante (Bk≠0) avec un gradient de température (ρF≠ρI).
La vitesse est passée par les arguments UX, UY et UZ. Les conditions aux limites de la vitesse sont COEFAX, COEFBX, ... Le flux de masse résultat est stocké dans les variables FLUMAS (faces internes) et FLUMAB (faces de bord). QDMX, QDMY et QDMZ sont des variables de travail qui serviront à stocker ρu, et COEFQA servira à stocker les Ã.
• Initialisation éventuelle du flux de masse
Si INIT vaut 1, le flux de masse est remis à zéro. Sinon, le
sous-programme rajoute aux variables FLUMAS et FLUMAB existantes le
flux de masse calculé.
• Remplissage des tableaux de travail
ρu est stocké dans QDM, et à dans COEFQA.
• Cas sans reconstruction
On calcule alors directement
FLUMAS=∑k=13[
α(ρI uk,I)+(1−α)(ρJ uk,J)]Sk
et
FLUMAB=∑k=13[
ρF Ak + BkρFuk,I]Sk
• Cas avec reconstruction
On répète trois fois de suite les opérations suivantes, pour k=1, 2 et 3
:
- Appel de GRDCEL
pour le calcul de grad (ρ uk).
- Mise à jour du flux de masse
FLUMAS=FLUMAS + [
α(ρI uk,I)+(1−α)(ρJ uk,J)
+1/2[grad (ρ uk)I+grad (ρ uk)J]
.OF]Sk
et
FLUMAB=FLUMAB+[
ρF Ak + BkρFuk,I
+Bkgrad (ρ uk)I.II′]Sk
• Annulation du flux de masse au bord
Quand le sous-programme a été appelé avec la valeur IFLMB0=1
(c’est-à-dire quand il est réellement appelé pour calculer un flux de
masse, et pas pour calculer le terme en div(ρR) par exemple), le flux
de masse au bord FLUMAB est forcé à 0, pour les faces de paroi et de
symétrie (identifiées par ISYMPA=0).