Dans ce sous-programme, appelé par codits et turbke , sont calculées les contributions au bilan explicite des termes convectifs et diffusifs reconstruits (sur maillages non orthogonaux et si l’utilisateur le désire) du second membre d’une équation de convection/diffusion pour un scalaire a. Ces termes s’écrivent de facon générique 1:
avec ρ, u, β et a variables connues à l’instant tn.
En s’inspirant des notations adoptées dans le sous-programme navsto , le bilan explicite correspondant à l’intégration sur une cellule Ωi de la partie convective −div( (ρ u)n a) de Bβ peut s’écrire sous forme d’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, s’ils existent. On a donc :
en posant :
où a f,ij et af bik représentent les valeurs aux faces
internes et de bord respectivement.
On rappelle :
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).
La valeur du flux convectif F ij dépend du type de schéma numérique choisi. Il en existe trois dans ce sous-programme :
| • un schéma décentré amont d’ordre 1 (upwind) : | |||||||||||||||||||||||||||||||||||||
| F ij((ρ u)n,a)= F ij amont((ρ u)n,a) | |||||||||||||||||||||||||||||||||||||
| où : | a f,ij=
| ||||||||||||||||||||||||||||||||||||
| • un schéma centré pur : | |||||||||||||||||||||||||||||||||||||
| F ij((ρ u)n,a)= F ij centré((ρ u)n,a) | |||||||||||||||||||||||||||||||||||||
| avec : | a f,ij= αij aI′+(1−αij) aJ′ | ||||||||||||||||||||||||||||||||||||
| • un schéma décentré amont d’ordre 2, SOLU 2(Second Order Linear Upwind) : | |||||||||||||||||||||||||||||||||||||
| F ij((ρ u)n,a)= F ij SOLU((ρ u)n,a) | |||||||||||||||||||||||||||||||||||||
| avec : | a f,ij=
| ||||||||||||||||||||||||||||||||||||
La valeur de F bik est calculée avec :
a bik est la valeur au bord donnée directement par les
conditions aux limites.
En schéma centré, on écrit en réalité (égalité conservant l’ordre 1 en espace sur a) :
On utilise le facteur 1/2 pour des raisons purement de stabilité numérique.
Un test de pente qui peut introduire des non linéarités dans l’opérateur de convection permet de basculer du schéma centré ou S.O.L.U. (d’ordre deux en maillage orthogonal) vers le schéma décentré amont d’ordre un (sans blending). De plus, en mode standard, on utilise en tout point une valeur de a f,ij issue d’une moyenne barycentrique entre la valeur décentrée amont et la valeur centrée (blending), suivant le souhait de l’utilisateur (variable BLENCV dans le sous-programme usini1 ).
De même, la partie diffusive peut s’écrire (au signe près) :
avec :
et :
en conservant les notations précédentes et avec S ij norme du vecteur
S ij, S bik norme du vecteur S bik,
a bik valeur au bord donnée directement par les conditions aux limites.
On rappelle le rôle des variables suivantes, intervenant lors de différents tests :
• IRCFLP, issu du tableau IRCFLU ; indique pour la
variable considérée si on reconstruit les flux de convection et de diffusion
= 0 : on ne reconstruit pas
= 1 : on reconstruit
• ICONVP, issu du tableau ICONV ; indique si on convecte
ou non la variable.
= 0 : on ne convecte pas
= 1 : on convecte
• IDIFFP, issu du tableau IDIFF ; indique si on prend en
compte la diffusion ou non pour la variable.
= 0 : on ne diffuse pas
= 1 : on diffuse
• IUPWIN indique localement dans bilsc2
(pour éviter
des calculs inutiles) si, pour la variable considérée à convecter, on choisit le schéma upwind pur (décentré amont) ou non.
= 0 : on n’est pas en upwind pur
= 1 : on travaille en upwind pur
• ISCHCP, issu du tableau ISCHCV ; indique, pour la variable considérée à
convecter, quel type de schéma convectif d’ordre deux sur maillage orthogonal
on choisit ( n’est utile que si BLENCP>0 ).
= 0 : on travaille avec le schéma S.O.L.U. (Second Order Linear
Upwind ) dit décentré amont d’ordre deux
= 1 : on travaille en centré
Dans ces deux cas, le coefficient de blending BLENCP est à fournir dans usini1
.
• BLENCP, issu du tableau BLENCV ; indique le pourcentage
de schéma convectif centré ou S.O.L.U. que l’on veut prendre en compte. Ce coefficient de pondération est compris entre 0 et 1.
• ISSTPP, issu du tableau ISSTPC ; indique si l’on veut
éviter le test de pente qui fait basculer le schéma convectif du second ordre
vers le schéma upwind si le test est positif.
= 0 : on travaille avec le test de pente systématiquement
= 1 : on travaille sans test de pente
Le gradient de la variable a traitée intervient dans le calcul du bilan
explicite. Un appel au sous-programme grdcel
qui permet de calculer cette quantité et
de la stocker dans le tableau (DPDX, DPDY, DPDZ) est donc
réalisé lorsque le calcul de ce gradient est nécessaire, c’est-à-dire :
• si l’option convection est activée avec un
schéma non upwind pur (ICONVP ≠ 0 et IUPWIN = 0) et,
si on veut reconstruire les flux (IRCFLP = 1),
ou si on veut utiliser le schéma S.O.L.U.
(ISCHCP = 0),
ou si on se sert du test de pente (ISSTPP = 0),
ou :
• s’il y a de la diffusion et si on veut
reconstruire les flux (IDIFFP ≠ 0 et IRCFLP = 1).
Sinon, le tableau (DPDX, DPDY, DPDZ) est mis à zéro.
On désigne par G c,i amont le gradient décentré de la
variable a, pour la
cellule Ωi. Il est stocké dans le tableau (DPDXA, DPDYA, DPDZA).
On définit également les scalaires a ij amont
et a bik amont par :
Après initialisation à
zéro, G c,i amont est calculé uniquement lorsque
l’utilisateur désire travailler en convection avec une méthode faisant
intervenir un schéma centré ou S.O.L.U. et en effectuant un test de pente.
• Pour chaque cellule Ωi, sont calculées les deux
quantités aIF (variable PIF) et aJF (variable PJF), valeurs à la face définies par :
Suivant le signe sijn du flux de masse (ρ u) ijn. S ij,
on affecte aIF ou aJF à la valeur a ij amont de
l’expression ∑j∈ Vois(i)a ij amont S ij.
• Quant aux termes de bord, ils sont calculés classiquement, en conservant les
notations adoptées dans les autres sous-programmes, comme suit :
Dans les tableaux (COEFAP, COEFBP) sont stockés (A b,ik,
B b,ik)k∈ γb(i), dans les tableaux (DIIPBX, DIIPBY,
DIIPBZ) le vecteur II′, dans SURFBO les surfaces
(S bik)k∈ γb(i).
Les contributions du bilan explicite [−div( (ρ u)na)
+div( β grad a)] sont calculées et ajoutées au tableau
second membre SMBR, qui a déjà été initialisé avant l’appel à
BILSC2
(par les
termes sources explicites par exemple, etc...).
La variable FLUX regroupe les parties convective et
diffusive des flux numériques. Elle se calcule classiquement sur les faces purement internes dans un
premier temps, puis sur les faces de bord. Les indices i et j sont
représentés respectivement par II et JJ.
Pour une prise en compte (lorsque nécessaire) aisée du signe sijn du
flux de masse (ρu) ijn. S ij, on utilise les relations
suivantes :
Pour un réel quelconque b, on a :
Dans ce sous-programme, b représente le flux de masse FLUMAS(IFAC)
pour une face interne IFAC (respectivement FLUMAB(IFAC) pour une
face IFAC de bord) ; b+ est stocké dans FLUI et b− dans FLUJ.
pour une face purement interne ij (IFAC)
Il s’agit de calculer :
La quantité sur laquelle porte la somme correspond à :
quelque soit le schéma convectif choisi, celui-ci impactant seulement les
affectations des quantités PIF (valeur à la face de a utilisée
quand b est positif) et PJF (valeur à la face de a utilisée
quand b est négatif). PIP représente aI′, PJP aJ′ et VISCF(IFAC)
β ij S ij/I′J′ .
La partie diffusive étant traitée de façon identique (soit en ne
reconstruisant pas, soit en reconstruisant), seul le schéma numérique
relatif à la convection diffère.
pour une face de bord ik (IFAC)
On s’occupe alors des termes :
avec :
Les coefficients ( A b,ik, B b,ik )k∈ γb(i) ( ( A b,ikdiff, B b,ikdiff )k∈
γb(i) ) traduisent les conditions aux limites associées à
a (respectivement au flux de diffusion 3 de a).
De même, la quantité sur laquelle porte la somme correspond à :
où PFAC représente a b1ik, PIP aI′, PFACD a bik et VISCB(IFAC)
β bik S bik/I′F .
Ce traitement est commun à tous les schémas, car les valeurs aux bords ne
sont fonction que des conditions aux limites et F bik est très
simplifié (upwind)4.
Il reste donc à déterminer, lorsque l’option convection est activée
(ICONVP = 1), les valeurs à affecter aux variables PIF et à
PJF, pour toute face IFAC interne commune à la cellule
Ωi et à la cellule Ωj.
Dans ce cas, aucune reconstruction n’est nécessaire puisque seules les valeurs
PVAR(II) et PVAR(JJ) au centre des cellules interviennent.
La variable INFAC comptabilise le nombre de calculs effectués en upwind
pur, pour impression dans le listing. Pour obtenir le flux numérique global FLUX (convectif + diffusif)
associé, on effectue les opérations suivantes :
• calcul des vecteurs II’ et JJ’,
• calcul du gradient facette (DPXF, DPYF, DPZF)
demi-somme des gradients volumiques G c,i et G c,j,
• calcul des valeurs éventuellement reconstruites aI′ et aJ′ (variables PIP et PJP respectivement) données par :
• calcul des quantités FLUI et FLUJ,
• calcul du flux FLUX via (??).
L’assemblage dans le second membre SMBR suit, en tenant compte de
l’expression (??)5 .
Les deux schémas possibles d’ordre deux sur maillage orthogonal sont ici soit le centré, soit le second ordre
décentré amont (S.O.L.U.).
Dans les deux cas, on effectue les opérations suivantes :
• calcul des vecteurs II’, tableau (DIIPFX, DIIPFY, DIIPFZ) et JJ’, tableau (DJJPFX, DJJPFY, DJJPFZ)
• calcul du gradient facette (DPXF, DPYF, DPZF)
demi-somme des gradients volumiques G c,i et G c,j,
• calcul des valeurs éventuellement (si IRCFLP = 1) reconstruites aI′ et aJ′ (variables PIP et PJP respectivement) données par :
• calcul des quantités FLUI et FLUJ.
sans test de pente (ISSTPP = 1)
en centré (ISCHCP = 1)
Les valeurs des variables PIF et PJF sont égales et calculées
à l’aide du coefficient de pondération αij comme
suit :
en S.O.L.U. (ISCHCP = 0)
Après avoir déterminé les vecteurs IF et JF, les valeurs
des variables PIF et PJF sont calculées comme suit :
On reconstruit systématiquement PIF et PJF afin d’éviter de
travailler en upwind pur, i.e. cette formule est appliquée même
lorsque l’utilisateur ne veut pas reconstruire (IRCFLP = 0).
avec test de pente (ISSTPP
= 0)
La démarche est identique à celle décrite dans le paragraphe
précédent, avec en plus, un calcul de test de pente qui fait rebasculer
localement (mais systématiquement) sous
certains critères le schéma centré ou S.O.L.U. choisi en schéma upwind pur.
calcul du test de pente
L’égalité (??) peut s’écrire, sur une cellule
Ωi purement interne, avec
sijn = sgn [(ρu) ijn . S ij] :
Sur une cellule Ωi de voisins (Ωj)j∈ Vois(i), le test de pente classique
consiste à repérer les non monotonies de la variable a en étudiant le signe du
produit scalaire des deux gradients volumiques G c,i
et G c,j. S’il est négatif, on rebascule en schéma upwind, s’il est
positif, on travaille en schéma centré ou S.O.L.U..
Une autre démarche garantissant également la monotonie de la solution est
d’appliquer ce critère aux gradients décentrés
G c,k amont ou à leur projection sur la normale à la
face ( G c,k amont . S kl ).
On étudie alors le signe du produit
G c,i amont . G c,j amont ou du produit
( G c,i amont . S ij ) . ( G c,j amont . S ij ) .
Le test de pente mis en œuvre est basé sur la première quantité,
G c,i amont . G c,j amont
(la seconde ayant été abandonnée, son champ d’action paraissant moins large). Le
choix d’un test basé sur
G c,i amont . G c,j amont
est tiré du raisonnement monodimensionnel suivant6 :
Soit p une fonction polynômiale de degré deux en x valant pI−1, pI,
pI+1 aux points I−1, I, I+1 d’abscisses respectives xI−1, xI et xI+1. Pour simplifier, on
supposera I en l’origine O ( xI = 0 ), le pas de discrétisation h constant,
soit xI+1 = − xI−1 = h . De plus, on suppose que la
vitesse est orientée du point I vers le point I+1, i.e. sijn =
1. C’est pourquoi on considère les points I−1, I et I+1 pour la face
ij qui se trouve ente I et I+1.
Le signe du produit p′(xI−1) . p′(xI+1) traduit la monotonie de la
fonction p. Si ce produit est positif, la fonction est monotone et on
travaille en centré ou en S.O.L.U., sinon, on repasse en schéma upwind. En identifiant les
coefficients du polynôme à l’aide des égalités p (xI−1) = pI−1
, p (xI) = pI , p (xI+1) = pI+1 , on obtient :
soit :
Or :
♣ pI+1 − pI/h représente la
dérivée décentrée au point I+1, directement accessible par les
valeurs de p dans les cellules voisines de la face ij,
♣ pI+1 − pI−1/2h la valeur de
la dérivée
centrée (au sens volumes finis) au point I, soit G c,i,
♣ pI − pI−1/h la valeur de la dérivée
décentrée (au sens volumes finis) au point I, soit G c,i amont
par définition.
Le test de pente relatif à l’expression p′(xI−1) . p′(xI+1) se
traduit donc par l’étude du signe de l’expression TP1d :
En raisonnant de façon analogue, une extension possible en dimension
supérieure consiste à remplacer les valeurs G c,k et G c,k amont par
( G c,k . S kl ) et
( G c,k amont . S kl ) respectivement. Ce qui
conduit à la formule TP3d+ :
ceci pour (ρ u) ijn. S ij > 0.
De même, on peut déduire une formulation TP3d−
associée à (ρ u) ijn. S ij < 0, définie par
:
Sont introduites les variables TESTI, TESTJ et TESTIJ qui valent :
La quantité TESQCK correspondant à TP3d, est
calculée dynamiquement, en fonction du signe du flux de masse sijn, afin de décentrer correctement.
alors :
⋄ si (ρ u) ijn. S ij > 0 et
si (G c,i . S ij)2 − (G c,i amont . S ij − a J−
a I/I′J′ S ij)TESQCK2 < 0 ou (G c,i amont . G c,j amont) < 0,
ou :
si (ρ u) ijn. S ij < 0 et
si (G c,j . S ij)2 −
(G c,j amont . S ij − a J−
a I/I′J′ S ij)TESQCK2 < 0 ou
(G c,i amont . G c,j amont) < 0,
alors on passe en upwind pur :
et INFAC est incrémenté.
⋄ sinon :
on affecte les valeurs en centré ou S.O.L.U. comme précédemment :
en centré (ISCHCP = 1)
Les valeurs des variables PIF et PJF sont égales et calculées
à l’aide du coefficient de pondération αij :
en S.O.L.U. (ISCHCP = 0)
Après avoir déterminé les vecteurs IF et JF, les valeurs
des variables PIF et PJF sont calculées comme suit :
On reconstruit systématiquement PIF et PJF afin
d’éviter de travailler en upwind pur, i.e. cette formule est appliquée même
lorsque l’utilisateur ne veut pas reconstruire (IRCFLP = 0).
Que l’on active le test de pente ou pas, lorsque l’option schéma centré ou
S.O.L.U. est activée, un coefficient de blending (BLENCP), fourni par l’utilisateur et compris
entre 0 et 1, permet de panacher si on le souhaite, le schéma considéré et le schéma décentré amont (upwind pur) suivant la formule :
• calcul des quantités FLUI et FLUJ,
• calcul du flux FLUX via (??).
L’assemblage dans le second membre SMBR suit, en tenant compte de
l’expression (??)7.
Pour plus d’information sur les schémas convectifs et le test de pente dans Code Saturne(version 1.1), on pourra se reporter au rapport EDF HI-83/04/020 (F. Archambeau, 2004).
• Schéma de convection
Schéma décentré amont (upwind pur)
Comme tout schéma d’ordre un, il est robuste, mais diffuse beaucoup numériquement.
Schéma centré ou S.O.L.U.
Cette classe de schémas peut générer des oscillations numériques, qui
peuvent provoquer l’explosion du calcul. Elle peut également générer
des dépassements non physiques sur des grandeurs scalaires physiques.
Compte-tenu de ces défauts et du non respect du principe du maximum, d’autres
schémas sont en cours de test et
d’intégration afin d’augmenter la qualité du
schéma proposé à l’utilisateur.
• Schéma de diffusion
La formule :
n’est d’ordre 2 que pour αij = 1/2.
Une correction à apporter pourrait être d’écrire :
avec un limiteur de gradient et un calcul de βij ne dégradant pas l’ordre.
• Réécriture informatique
Afin d’améliorer le temps C.P.U., un effort peut être fait au niveau de l’écriture des boucles. Notamment, il existe un test IF dans une boucle sur la variable IFAC qu’il faudra regarder.
• Calcul du gradient intervenant lors de la reconstruction des flux diffusifs
Pourquoi prend-on l’expression 1/2 ( G c,i + G c,j ) au lieu de G c,k , pour k=i ou pour k=j dans les valeurs reconstruites aI′ ou aJ′ de (??) et (??)?