Le but de ce sous-programme, appelé par codits
et covofi
, est de construire la
matrice de convection/diffusion, incluant les contributions adéquates des termes sources,
intervenant dans le membre de gauche d’équations discrétisées telles que
celle de la
quantité de mouvement, d’une équation de convection diffusion d’un scalaire
ou de modèle de turbulence.
Le type d’équation considérée est, pour la variable scalaire a :
La matrice ne s’applique qu’aux termes non reconstruits, les autres étant pris en compte au second membre et
traités dans le sous-programme bilsc2
. La partie
convective, lorsqu’elle existe, est issue du schéma upwind pur, quelque soit
le type de schéma convectif choisi par l’utilisateur. En effet, c’est, à
l’heure actuelle, la seule façon d’obtenir un opérateur linéaire à
diagonale dominante.
La matrice est donc associée à EMscal, opérateur
agissant sur un scalaire a (inspiré de celui vectoriel EM
défini dans navsto
) d’expression actuelle, pour tout cellule Ωi de
centre I :
avec :
• fsimp le coefficient issu du terme instationnaire
ρ |Ωi|/Δ t, s’il y a lieu, et de
l’implicitation de certains termes sources (contribution découlant de la prise
en compte de la
variation ∂ ρ /∂ t de
la masse volumique ρ au cours du temps, diagonale du tenseur de pertes de
charges par exemple...): cette initialisation est en fait effectuée en amont
de ce sous-programme,
• F ij amont le flux numérique convectif scalaire
décentré amont calculé à la face interne ij de la cellule Ωi,
• F bik amont le flux numérique convectif scalaire
décentré amont associé calculé à la face de bord ik de la cellule Ωi
jouxtant le bord du domaine Ω,
• D ij NRec le flux numérique diffusif scalaire non
reconstruit associé calculé à la face interne ij de la cellule Ωi,
• D bik NRec le flux numérique diffusif scalaire non
reconstruit associé calculé à la face de bord ik de la cellule Ωi jouxtant le bord du domaine Ω,
• Vois(i) représente toujours l’ensemble des cellules
voisines de Ωi et γb(i) l’ensemble des faces de
bord de Ωi.
Du fait de la résolution en incréments, a est un incrément et ses
conditions aux limites associées sont donc de type Dirichlet homogène ou de
type
Neumann homogène.
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).
L’opérateur EMscal s’écrit, pour tout I centre de cellule :
où
a f,ij = a I ou a J
selon le signe de (ρ u) ijn.S ij (schéma
convectif upwind systématique),
et avec I′J′, mesure algébrique, orientée comme la
normale sortante à la face, i.e. allant de I vers J pour la cellule
Ωi de centre I. On la
notera I′J′ I lorsqu’on aura besoin d’expliciter
clairement l’orientation.
af bik = aI ou
a bik selon le signe de
(ρ u) bikn. S bik (schéma
upwind systématique)
et a bik valeur au bord est donnée directement par les conditions
aux limites (valeur non reconstruite). I′F, mesure algébrique, orientée relativement à la
normale sortante à la face, i.e. allant de I vers l’extérieur du domaine.
En général, sauf cas pathologiques, les mesures algébriques
I′J′ et I′F
sont positives et correspondent aux distances I′J′ et I′F. On se reportera
au paragraphe Points à traiter pour plus de détails.
Soit EM scal la matrice associée ; sa taille est a priori de
NCEL * NCEL, mais compte-tenu de la nature de la structure de
données, seuls deux tableaux DA(NCEL) contenant les valeurs
diagonales et XA(NFAC,*) les contributions des termes extra-diagonaux sont nécessaires, avec NCEL nombre de
cellules du maillage considéré et NFAC nombre de faces internes associé.
Du fait des simplifications effectuées sur la matrice (non reconstruction des
termes), les composantes extradiagonales de la ligne I ne sont différentes de zéro que pour
les indices J des cellules voisines de I. On peut donc stocker toutes les
contributions non nulles de la matrice dans deux tableaux DA(NCEL) et XA(NFAC,2) :
• DA(I) est le coefficient de la colonne I dans la ligne I,
• si IFAC est une face qui sépare les cellules Ωi
et Ωj, orientée de I vers J, alors :
XA(IFAC,1) est le coefficient de la colonne J dans la ligne I et
XA(IFAC,2) est le coefficient de la colonne I dans la ligne J.
Lorsque la matrice est symétrique, i.e. lorsqu’il n’y a pas de convection (ICONVP = 0) et que seule la
diffusion est à prendre en compte, alors XA(IFAC,1) = XA(IFAC,2)
et on réduit XA à XA(NFAC,1).
Soit m ijn ( m bikn ) la valeur de (ρ
u) ijn.S ij (respectivement (ρ
u) bikn. S bik).
Alors :
contribution volumique : fs imp aI
contribution d’une face purement interne ij
L’expression
s’écrit :
Ici, I′J′ = I′J′ I.
contribution d’une face de bord ik
De même :
avec :
a n’étant associée qu’à des conditions aux limites de type Dirichlet
homogène ou de type Neumann homogène.
L’indicateur de symétrie ISYM de la matrice considérée est affecté comme suit :
• ISYM = 1 , si la matrice est symétrique ;
on travaille en diffusion pure , ICONVP = 0 et IDIFFP = 1,
• ISYM = 2 , si la matrice est non symétrique ;
on travaille soit en convection pure ( ICONVP = 1, IDIFFP = 0 ), soit en
convection/diffusion ( ICONVP = 1, IDIFFP = 1 ).
Les termes diagonaux de la matrice sont stockés dans le tableau
DA(NCEL). Ceux extra-diagonaux le sont dans XA(NFAC,1) si la
matrice est symétrique, dans XA(NFAC,2) sinon.
Le tableau DA est initialisé à zéro pour un calcul avec
ISTATP = 0 (en fait, ceci ne concerne que les calculs relatifs à la
pression). Sinon, on lui affecte la valeur ROVSDT comprenant la partie instationnaire, la contribution du terme continu en − a div(ρ u)n
et la partie diagonale des termes sources implicités. Le tableau
XA(NFAC,*) est initialisé à zéro.
Ils ne se calculent que pour des faces purement internes IFAC, les faces de bord ne contribuant qu’à la diagonale.
Pour chaque face interne IFAC, les contributions extradiagonales relatives
au terme aI et à son voisin associé aJ sont calculées dans
XA(IFAC,1) et XA(IFAC,2) respectivement (pour une face orientée
de I vers J).
On a les relations suivantes :
avec FLUMAS(IFAC) correspondant à m ijn, FLUI à 1/2 ( m ijn − |
m ijn| ), VISCF(IFAC) à β ij
S ij/I′J′ .
XA(IFAC,1) représente le facteur de aJ dans la
dernière expression de (??).
FLUJ correspond à −1/2 ( m ijn + |
m ijn| ). En effet, XA(IFAC,2) est le facteur de aI dans l’expression équivalente de
la dernière ligne de (??), mais écrite en J.
Ce qui donne :
Le terme recherché est donc :
1/2( m jin − | m jin| )−
β jiS ji/J′I′ .
Or :
m jin = − m ijn (S ji = −
S ij et (ρ u) jin = (ρ u) ijn
), avec J′I′ mesure algébrique, orientée relativement à la
normale sortante à la face, i.e. allant de J vers I. On la note
J′I′ J.
On a la relation :
d’où la deuxième égalité dans (??).
Pour chaque face interne IFAC, la contribution extradiagonale relative au
terme aI est calculée dans
XA(IFAC,1) par la relation suivante :
avec VISCF(IFAC) à β ij
S ij/I′J′ .
Pour chaque face interne ij ( IFAC ) séparant les cellules Ωi
de centre I et Ωj de centre J, DA(II) est la quantité en facteur de aI dans la
dernière expression de (??), soit :
De même, pour DA(JJ), on a :
d’après l’expression de (??) et l’égalité (??).
L’implantation dans Code Saturneassociée est la suivante :
pour toute face IFAC d’éléments voisins II =
IFACEL(1,IFAC) et JJ = IFACEL(2,IFAC),
on ajoute à
DA(II) la contribution croisée −XA(IFAC,2) (cf.
(??)) et
à
DA(JJ) la contribution −XA(IFAC,1) (cf.
(??)).
Elles interviennent juste dans le tableau DA, compte-tenu de leur
écriture et définition. Elles se calculent via la dernière
expression de (??). Pour chaque face IFAC, de l’élément
de centre I, jouxtant le bord, on s’intéresse à :
avec :
soit :
ce qui, pour le terme sur lequel porte la somme, se traduit par :
ICONVP * (− FLUJ + FLUI * COEFBP(IFAC) + IDIFFP *
VISCB(IFAC) * ( 1 − COEFBP(IFAC))
avec,
m bikn représenté par FLUMAB(IFAC) ,
1/2 (
m bikn + | m bikn| ) par - FLUJ ,
1/2 ( m bikn −
|m bikn|)B b,ik par FLUI ,
B b,ik par COEFBP(IFAC), β bikS bik/I′F par VISCB(IFAC).
Lorsqu’il n’existe aucune condition à la limite de type Dirichlet et que ISTATP = 0 (c’est-à-dire pour la pression uniquement), on déplace le spectre de la matrice EM scal de EPSI (i.e. on multiplie chaque terme diagonal par (1 + EPSI) ) afin de la rendre inversible. EPSI est fixé en dur dans matrix à 10−7.
• Initialisation
Le tableau XA est initialisé à zéro lorsqu’on veut annuler la
contribution du terme en
ρ |Ωi|/Δ t, i.e. ISTATP = 0 . Ce qui ne permet donc pas la prise en
compte effective des parties diagonales des termes sources à impliciter,
décidée par l’utilisateur. Actuellement, ceci ne sert que pour la variable
pression et reste donc a priori correct, mais cette démarche est à
corriger dans l’absolu.
• Nettoyage
La contribution ICONVP FLUI, dans le calcul du terme
XA(IFAC,1) lorsque la matrice est symétrique est inutile, car
ICONVP = 0.
• Prise en compte du type de schéma de convection dans
EM scal
Actuellement, les contributions des flux convectifs non reconstruits sont
traitées par schéma décentré amont, quelque soit le schéma choisi par
l’utilisateur. Ceci peut être handicapant. Par exemple, même sur
maillage orthogonal, on est obligé de faire plusieurs sweeps pour obtenir une
vitesse prédite correcte. Un schéma centré sans test de pente peut
être implanté facilement, mais cette écriture pourrait, dans l’état
actuel des connaissances, entraîner des instabilités
numériques. Il serait souhaitable d’avoir d’autres schémas tout aussi
robustes, mais plus adaptés à certaines configurations.
• Maillage localement pathologique
Il peut arriver, notamment au bord, que les mesures algébriques,
I′J′ ou I′F soient négatives (en cas de maillages
non convexes par exemple). Ceci peut engendrer des problèmes plus ou moins
graves : perte de l’existence et l’unicité de la solution (l’opérateur associé n’ayant plus les bonnes propriétés de régularité
ou de coercivité), dégradation de la matrice (perte de la positivité) et donc résolution par solveur linéaire
associé non approprié (gradient conjugué par exemple).
Une impression permettant de signaler et de localiser le problème serait souhaitable.