BRGM MEFISTO STEF Manuel utilisateur

application de la méthode des éléments finis

à la simulation des transferts dans les eaux souterraines

notice et m o d e d'emploi du programme M E F I S T O - S T E F

î. S / R * .

BRGM

application de la méthode des éléments finis

à la simulation des transferts dans les eaux souterraines

notice et mode d'emploi du programme MEFISTO-STEF

M. Récan octobre 1986

86 SGN 194 EAU

BUREAU DE RECHERCHES GÉOLOGIQUES ET MINIÈRES

SERVICE GÉOLOGIQUE NATIONAL

Département Eau

B.P. 6009 - 45060 ORLÉANS CEDEX 2 - Tél.: 38.64.34.34

RESUME

Le présent rapport comprend la description, le m o d e d'emploi et un jeu d'essai du programme M E F I S T O - S T E F (STEF = simulation des Transferts par Eléments

Finis) dans la version transfert hydrodispersif (déplacement d'un soluté avec dispersion, transfert de chaleur par convection forcée...).

L'utilisation de ce code de calcul nécessite un minimum de connaissances sur les méthodes de calcul par éléments finis ; celles-ci pourront être acquises par une lecture approfondie de l'ouvrage de Dhatt et Touzot, qui fournit un programme de base à partir duquel M E F I S T O et M E F I S T O - S T E F ont été conçus. Tous les apports nouveaux par rapport à ce programme sont détaillés dans le présent document.

Le choix de la discrétisation et l'élimination des problèmes de stabilité, convergence... demande une expérience solide qui ne peut être acquise et conservée que par l'utilisation régulière de ces techniques. La mise en oeuvre de ce logiciel devra être réservée aux problèmes difficiles requiérant une précision élevée.

Le programme M E F I S T O a été mis au point à l'occasion d'études diverses sur contrat ; des compléments, la généralisation du code de calcul, et la préparation du m o d e d'emploi ont été réalisés sur fonds propres du Ministère de l'Industrie dans le cadre du programme Informatique Hydrogéologique du Département E A U .

Les principales études réalisées à l'aide de M E F I S T O sont :

- différents problèmes de stockage d'eau chaude dont celui du pilote de l'Ecole Normale Supérieure (Lyon Gerland),

- la simulation de la gazéification souterraine du charbon (rétrocombustion),

- les transferts de substance dissoute en aval d'un stockage souterrain de gaz naturel,

- pompage par un puits incomplet.

SOMMAIRE

Page

INTRODUCTION 1

1 - LA BIBLIOTHEQUE D'ELEMENTS CV2D 2

1.1 - L'équation de diffusion-convection 2

1.2 - Notations et conditions aux limites 4

1.3 - Formulation variationnelle 6

1.3.1 - Formulation en flux total 7

1.3.2 - Formulation en flux diffusif 8

1.4 - Discrétisation spatiale et formulation matricielle 10

1.5 - Calcul effectif des matrices élémentaires 12

1.5.1 - Transformation des dérivées premières 12

1.5.2 - Transformation des intégrales 13

1.5.3 - Intégration numérique 13

1.5.4 - Organisation du calcul 13

1.5.5 - Calcul de la matrice de rigidité 15

1.5.6 - Calcul de la matrice de masse 18

1.5.7 - Calcul des sollicitations réparties 19

1.5.8 - Calcul des sollicitations sur les frontières 19

1.6 - Utilisation de la bibliothèque C V 2 D 21

2 - NOTIONS S U R LES S C H E M A S U P W I N D , L A STABILITE ET L A

DISPERSION DES S C H E M A S N U M E R I Q U E S 23

2.1 - Les schémas "upwind" 23

2.1.1 - C a s monodimensionnel 23

2.1.2 - Extension à 2 ou 3 dimensions 26

Page

2.2 - Notions de stabilité, dispersion et dissipation des schémas numériques 27

2.2.1 - Stabilité (équation de la diffusion) 27

2.2.2 - Dissipation, dispersion 31

2.3 - Oscillations avec l'équation de la diffusion 35

2.4 - L e principe du m a x i m u m 38

2.4.1 - Définition 38

2.4.2 - Cas de l'équation de la diffusion 39

2.4.3 - Cas de l'équation de diffusion-convection monodimensionelle 43

2.5 - Choix du schéma de discrétisation 43

2.5.1 - Les causes des oscillations numériques 43

2.5.2 - L e nombre de Péelet numérique 44

2.5.3 - Les schémas upwind 44

2.5.4 - Diagonalisation de la matrice masse 45

2.5.5 - Choix des éléments 45

2.6 - Exemple 46

Références bibliographiques 48

INTRODUCTION

D e nombreux problèmes physiques sont modélisés par l'équation de diffusion convection. Mis à part le domaine de la mécanique des fluides qui fait intervenir l'équation de diffusion-convection avec des termes supplémentaires non linéaires (inertie), cette équation est à la base des problèmes de pollution ou de transferts thermiques dans les aquifères.

La méthode des éléments finis a connu un essor considérable ces vingt dernières années. Elle est à la base des codes de calcul utilisés en mécanique des solides et s'est avérée très efficace dans le traitement des problèmes basés sur l'équation de la diffusion. Les avantages de la méthode des éléments finis résident dans la possibilité de représenter facilement des domaines de forme géométrique compliquée et de traiter de façon consistente les conditions aux limites. Sa mise en oeuvre est cependant sensiblement plus compliquée que dans le cas des différences finies et elle coûte en général plus chère.

C e rapport présente la bibliothèque d'éléments CVZD basée sur la formulation de Galerkine. Couplée au programme M E F T S T O (1), cette bibliothèque permet de traiter un problème de diffusion convection. L a formulation utilisée et le m o d e d'emploi de la bibliothèque CVZD sont exposés dans le premier chapitre. Sa lecture suppose un minimum de connaissances de la méthode.

L'application de la méthode de Galerkine au traitement de l'équation de diffusion-convection fait apparaître des problèmes similaires à ceux rencontrés avec la méthode des différences finies lorsque les transferts convectifs sont importants.

Le deuxième chapitre présente quelques notions de base sur la stabilité et la dispersion des schémas à partir d'exemples simples et m e t en évidence les causes des oscillations qu'on observe parfois dans les solutions numériques.

U n exemple d'utilisation de la bibliothèque CVZD est proposé dans le dernier chapitre.

1 - LA BIBLIOTHEQUE D'ELEMENTS CV2D

1.1 - L'EQUATION DE DIFFUSION-CONVECTION

En adoptant la terminologie associée au traitement des problèmes de pollution, l'équation de diffusion-convection s'écrit :

dC * ••• • •

w — + V({J.C) = V

dt

D VC avec C : concentration to : porosité cinématique du milieu

U : vitesse de Darcy de l'écoulement

Q : terme source

D : tenseur de diffusion/dispersion.

(1.1)

C e tenseur est lui-même la s o m m e de deux tenseurs :

- D o

= do.n, d

0

étant le coefficient de diffusion moléculaire et

II le tenseur identité

- D a

: tenseur de dispersion supposé symétrique dont le premier vecteur propre a la direction du vecteur vitesse de l'écoulement et le second espace propre est constitué par le plan orthogonal au vecteur vitesse. E n géométrie plane et dans un repère lié à l'écoulement, D a

s'écrit :

D = I

a

a

L 0

0

a

L

\u\

où a L est la dispersivité longitudinale (dans le sens de l'écoulement) a-p est la dispersivité transversale (dans le plan perpendiculaire à l'écoulement)

IUI est le module d u vecteur vitesse.

Si U x

et U y sont les composantes du vecteur vitesse, dans le repère du problème, les composantes du tenseur de diffusion/dispersion dans ce repère s'écrivent :

D = xx

a, U

2

+a

T

U

2 yy \u\ °

(1.2)

(a, -a_)î/ U

D = D = —

X

—^

*y y* \u\

It/I = ( t /

2

+

1 j

L'équation (1.1) permet, à condition de réinterpréter les coefficients, de traiter les problèmes de transferts thermiques dans un aquifère. L'équation régissant les transferts thermiques s'écrit :

E

?¿/

T )

À

VT

(1.3)

Y f

dt

avec T : température

Q-p : terme source

U : vitesse de Darcy de l'écoulement

Y a

: capacité calorifique du milieu (tenant compte de la présence de la phase fluide)

Yf : capacité calorifique de la phase fluide

X : tenseur de conductivité équivalente qui regroupe la conductivité isotrope du milieu poreux (solide + fluide) en l'absence d'écoulement et un terme de dispersion lié à l'hétérogénéité de la vitesse, fonction linéaire de cette vitesse. Dans les axes longitudinaux et transversaux liés à la vitesse, les deux composantes diagonales du tenseur X sont :

(1.4) où B L et 6 T

s o n t

respectivement les dispersivités thermiques longitudinales et transversales.

L'équation (1.3) est formellement identique à l'équation (1.1) à condition d'identifier w et Y a

/Yf> ¿o

e t

^o/ïf»

a

L

e t a

T> &L

e t

BT> Q et

1.2 - NOTATIONS ET CONDITIONS AUX LIMITES

La bibliothèque d'éléments C V 2 D est destinée à traiter un problème de diffusion-convection en géométrie plane ou axisymétrique. Dans le développement qui suit, nous supposerons que le nombre de dimension de l'espace est n¿ ce qui permet de donner un aspect tout à fait général à la formulation.

Soit V un domaine ouvert de R

1 1

^ et S sa frontière. Soit x = {xj}, i = 1, 2,

..., n¿ un point de V et n = {nj} la normale unitaire extérieure à S. Nous supposons que la frontière S peut être décomposée de la façon suivante :

S = S U S

C <J

(1.5) cj) = s n s où S c

et Sq sont des sous-ensembles de S (dont l'intersection est vide et la réunion est

S). La barre superposée indique la fermeture de l'ensemble et 0est l'ensemble vide.

Nous utiliserons la convention de sommation implicite sur les indices répétés (par exemple si n = Z, alors qjnj = q^n^ + q2 n

2)« U n e virgule est utilisée pour indiquer une dérivtion partielle (par exemple qi i = oqj/Sxj).

Considérons une triangulation du domaine V en éléments V e

, e = 1, 2, ..., rig où n e

est le nombre d'éléments. L a frontière de chaque élément est notée S e

et nous supposons

_

U

V =

u

V e

(1.6)

S C

u

S e

U e

représente les intérieurs des éléments (hors frontières) et l'ensemble Sjnt = U S e

- S représente les frontières intérieures de la triangulation.

Se

L'équation (1.1) fait apparaître 3 types de flux

- le flux convectif q

c

= U . C (1.7)

I I

- le flux diffusif 9? = -

D a

C

>j Í

1

-

8

)

- le flux total <?, = q* + q'f (1.9)

L'équation (1.1) peut, compte tenu de ces notations, s'écrire sous la forme :

q . = f= Q -

w

dC/dt

(1.10)

Les notations suivantes seront utiles dans la suite

V =V..n. q = q..n q c

= q c

. n. q d

= q d

n. (1.11)

C o m p t e tenu de (1.10), nous ne distinguerons pas dans un premier temps le cas stationnaire du cas instationnaire (le cas instationnaire requiert cependant l'introduction de conditions initiales).

Le but du problème est de trouver la fonction C satisfaisant l'équation

(1.10) ainsi que certaines conditions aux limites.

Nous imposons C = g sur S (1.12)

c

où g est une fonction donnée définie sur S c

. L a portion S c

de la frontière S est une frontière à concentration imposée. Il reste à définir une condition sur Sq. Nous envisagerons Z types de conditions

- q

n

= l sur S (1.13)

-q^hsurS (1-1*) où f et h sont des fonction données définies sur Sq.

La première condition porte sur le flux total pénétrant dans le domaine par la frontière Sq et la seconde sur le flux diffusif pénétrant dans V .

R e m a r q u e : Bien que ce ne soit pas toujours précisé, la plupart des formulations développées pour l'équation de diffusion-convection satisfont la condition

(1.14) portant sur le flux diffusif.

1.3 - FORMULATION VARIATIONNELLE

La méthode des résidus pondérés consiste à rechercher des fonctions d'interpolations C qui annulent la forme intégrale issue de (1.10).

V iq -f).

Les fonction d'interpolation sont supposées satisfaire

C = g sur S

et les fonctions de pondération

W = 0 sur S

(1.16)

(ce qui est en pratique réalisé en modifiant le système d'équation résultant de la discrétisation en éléments finis).

L'étape suivante consiste à diminuer l'ordre maximal des dérivés qui interviennent dans (1.15) (afin de diminuer les conditions de dérivabilité des fonctions d'interpolation). Cette opération est réalisée par le théorème d'Ostrogradski qui

énonce :

= - f V^

Etant donné que seul le flux diffusif fait intervenir les dérivées à l'ordre 1 des fonctions de pondération, on peut, pour diminuer l'ordre maximal de dérivation, appliquer le théorème d'Ostrogradski

- soit au terme en flux dif fusif seulement

- soit au terme en flux total.

Ces deux possibilités conduisent à des hypothèses de nature différentes sur la continuité des flux à l'intérieur du domaine et sur les conditions aux limites qui sont prise en compte de façon naturelle.

suivante :

Avant d'expliciter les deux formulations, nous introduisons la notation

Soit x un point de S m t, c'est-à-dire situé sur une frontière interne de la triangulation. O n appelle (de façon arbitraire) un côté de S^ le côté positif et l'autre

le côté négatif. Soient n

+

et n~ les normales unitaires à S m t

au point x qui sont dirigées vers les directions positives et négatives. Il est évident que n

+

= - n~.

Soient qj

+

et qj~ les valeurs de qj obtenues en approchant du point x par les côtés positifs et négatifs. L e saut de q n

au point x est défini par :

- \q.

+

— q. n . = q.

+

n.

+

+ q. n.

(1.18)

II est clair que le saut [q n

] est invariant si l'on intervertit les côtés positifs et négatifs.

1.3.1 - FORMULATION EN FLUX TOTAL

Pour définir l'équation variationnelle, nous supposons que toutes les fonctions (et en particulier les fonctions C et W ) sont de classe C

0 0

(i.e. indéfiniment continues et dérivables).

E n appliquant le théorème d'Ostrogradski à l'équation (1.5), on obtient :

- I q

i

W , . d V = f W d V -

q i t

(1.19)

Décomposons l'intégrale sur la frontière du domaine finalement

qn WdS= q WdS + q WdS

S

l l

¡S

n

¡S

n c q

Puisque W est nulle sur S c

et compte tenu de (1.13), ce terme s'écrit

q.n.WdS = - W. l. dS

S

l

' ¡S

Q

et l'équation variationnelle équivalente à (1.15) est :

- qW dV = fWdV+ f W.ldS

(

->l>

* V ¡S

1

Dans la méthode des éléments finis, les intégrales portant sur le domaine sont remplacées par la s o m m e des intégrales sur chaque élément de la triangulation.

Si les fonctions W et C sont de classe C

0 0

, l'équation variationnelle (1.20) est équivalente à :

C o m m e les fonctions de pondération sont arbitraires, les équations d'Euler-Lagrange associées sont bien (1.10) et (1.13).

Si W et C sont seulement de classe C ° (et plus exactement telles que leurs dérivées soient des continues sur les frontières internes de la triangulation, i.e. sur

Sint) l'équation variationnelle (1.20) est équivalente à

s. • ~ . .

h int

L'équation (1.22) montre que, dans la classe des fonctions C ° , les

équations d'Euler-Lagrange sont

- la restriction à l'intérieur des éléments de (1.10)

- l'équation (1.13)

- la condition de continuité du flux total à travers les frontières internes de la triangulation, c'est-à-dire :

[qj - 0 sur

int ( U 2 3 )

7.3.2 - FORMULATION EN FLUX DIFFUSIF

Dans cette formulation, le théorème d'Ostrogradski n'est appliqué qu'au terme de flux diffusif. E n ré-écrivant (1.15) sous la forme :

*. +q

c¡t

- f)WdV = 0

et en tenant compte de :

WdV=- q d

,WdV+ q d

Wn dS

1 l

h

l

on obtient :

q d

\V,idV+ q c

. W d V = WfdV- q

d n WdS

9 c

soit en tenant compte de (1.14) et puisque W est nulle sur S c

:

q d

W .dV + q c

I

S

Q

Si les fonctions W et C sont de classe C

0 0

, l'équation variationnelle (1.24) est équivalente à

W(q -f)dV-\ W(q d

+h)dS = 0 (1.25)

¡u

l

'

1

¡S,

n

n

et les équations d'Euler-Lagrange sont (1.10) et (1.14).

Si W et C ne sont que de classe C ° , l'équation (1.24) conduit à :

V I W(q

t

.-f)dV-\ W(q d

+h)dS-\ W[q d

] dS = 0 (1.26)

e

h int

Les équations d'Euler-Lagrange sont alors

- la restriction de l'équation (1.10) à l'intérieur des éléments

- l'équation (1.14)

- la condition de continuité du flux diffusif à travers les frontières internes de la triangulation, c'est-à-dire

= 0 sur S. (1.27)

int

II apparaît finalement que le choix d'une formulation variationnelle ((1.20) ou (1.24)) conduit, dans le cadre des fonctions de classe C ° (ce qui est le cas des

éléments usuels), à des conditions différentes.

La formulation en flux total permet de satisfaire naturellement les conditions aux limites en flux total imposées sur la frontière S et implique la continuité du flux total à travers les frontières des éléments. Dans la formulation en flux diffusif, les conditions sont équivalentes mais s'appliquent au flux diffusif.

10

Le choix d'une formulation est largement guidé par le type de conditions que l'on veut imposer sur la frontière du domaine d'étude. Dans le cadre des problèmes susceptibles d'être traités au Département E A U , les conditions aux limites porteront plutôt sur le flux diffusif. Par ailleurs, il est souvent nécessaire de simuler une condition d'exhaure sur une portion de la frontière. Généralement, cette frontière est placée, de façon un peu arbitraire, dans une région où les gradients sont faibles et loin de la zone qui présente le plus d'intérêt. Sur le plan numérique, il est souhaitable que la condition d'exhaure n'introduise pas un "bruit numérique" susceptible d'être propagé vers la zone d'intérêt. Dans les problèmes dominés par la convection, il apparaît que la condition en flux diffusif est plus satisfaisante à cet égard. En effet, si le terme de diffusion est faible par rapport au terme de convection, la condition de flux total nul sur une frontière correspondant à un exhaure donné :

0 = - q = -U C+ q

d

= -U C ce qui tend à imposer C = 0. C e type de condition tend à créer une couche limite au voisinage de la frontière et génère donc souvent des difficultés numériques (cf. Z.l).

Par contre, une condition en flux diffusif nul à l'exhaure implique (si D est isotrope)

8 C / 9 n = 0 ; ce qui apparaît cohérent avec l'hypothèse que l'exhaure se trouve dans une région où les gradients sont très faibles.

1.4 - DISCRETISATION SPATIALE ET FORMULATION MATRICIELLE

Nous utiliserons la formulation en flux diffusif associée à la méthode de la

Galerkine, c'est-à-dire que les fonctions de pondération sont égales aux fonctions d'interpolation de la solution sur le domaine.

Dans la méthode des éléments finis, les intégrales sur le domaine sont remplacées par la s o m m e des intégrales sur chaque élément. Par ailleurs, sur un

élément, les fonctions d'interpolation associées à un noeud n'appartenant pas à l'élément sont identiquement nulles. C'est ce qui explique l'aspect répétitif de la méthode qui a largement contribué à son succès. Il suffit alors de discrétiser l'équation variationnelle (1.24) sur un élément général.

Afin d'éviter une surcharge en indices, nous utiliserons l'opérateur de dérivation V et nous restreindrons l'espace à deux dimensions (qui constitue le champ d'application de la bibliothèque C V 2 D ) .

1 ]

Soit e un élément de la triangulation possédant n noeuds. La solution approchée sur cet élément s'écrit :

n

C=C N. {= Y C.N )

(

j j ¿— j J

Les coefficients C ; sont les valeurs de la solution aux noeuds et sont fonctions du temps pour un problème transitoire. Les fonctions d'interpolation N ; sur l'élément ne dépendent que des coordonnées spatiales des noeuds de l'élément.

L'équation variâtionnelle (1.24) s'écrit sur cet élément :

- | q

d

V N dV + \ .Vq

c

.N.dV = (Q - W BCloí) N. dV + N .hdS

¡V e l

¡V e l

¡v e

' ¡s

l

Q

C o m p t e tenu de la définition des flux ((1.7) et (1.8)) et de (1.28), l'équation précédente donne :

wN.N dV C + VN.DVN.dV+\ N.V(UN.)dV

e i

C. = QN dV +

(1.29)

S q

h.N dS

Cette équation peut s'écrire sous la forme matricielle suivante : avec {Ci} : valeurs nodales de la concentration

{ q } = d{Cj}/dt

[m] - \ e

w

^i^j } _ matrice de masse

[k] =

f

I

e v N

,

D v N

j

dV

matrice de diffusion/dispersion

) dV matrice de convection

{f v

} = r \ QN.dV vecteur des sollicitations réparties

{f s

} = j ' i ' vecteur des flux diffusifs imposés sur les frontières S n

Remarque : - E n géométrie plane (x, y) d V = dx.dy.l

dS = ds.l

- E n géométrie axisymétrique (r, z) dV = Et r dr dz dS = 2f r ds où s est l'abscisse curviligne le long de la frontière.

12

1.5 - CALCUL EFFECTIF DES MATRICES ELEMENTAIRES

La bibliothèque C V 2 D comporte à l'heure actuelle 5 types d'éléments isoparamétriques dont les fonctions d'interpolation sont de classe C ° . Le calcul des matrices et vecteurs élémentaires définis dans le paragraphe précédent est réalisé par intégration numérique sur l'élément de référence. Nous noterons Ç et ri les coordonnées spatiales sur l'élément de référence. Puisque les éléments sont isoparamétriques, les fonctions d'interpolation géométriques sont identiques aux fonctions d'interpolation de la variable sur l'élément. Si (Ç, n) est un point de l'élément de référence, le point de l'élément réel qui lui correspond par la transformation géométrique est (x, y) avec :

x — A M Ç , n)x.

(1.31)

y = NA(.n)y.

où (XJ, yj) sont les coordonnées des n noeuds de l'élément réel.

Pour ramener le calcul des matrices sur l'élément de référence, il faut transformer les dérivations qui interviennent dans la matrice de rigidité ainsi que les intégrales.

Í.5.1 - TRANSFORMATION DES DERIVEES PREMIERES

A partir de la règle de dérivation en chaîne qui permet de calculer les dérivées en Ç d'une fonction à partir de ses dérivées en x, on obtient : a

n

a*

dx

an

îl.

il J

a n

3

f — i

\±\-

[J] dy

a

Sx

a

~y

(1.32)

où [J] est la matrice jacobienne de la transformation géométrique.

Les termes de [J] sont obtenus en dérivant par rapport à Ç et X] la relation

(1.31). O n obtient par exemple :

Bx 9N.

— rr • •

Y

ai, dt x i

13

7.5.2 - TRANSFORMATION DES INTEGRALES

Le changement de variables défini par la relation (1.31) permet de passer de l'intégration d'une fonction f sur l'élément réel V e

à une intégration plus simple sur l'élément de référence V

1

"

V e

f(x).dx.dy= \

T

f{x{t,,n)).det{J).dl.dr\ (1.33) det (J) étant le déterminant de la matrice jacobienne [J].

7.5.3 - INTEGRATION NUMERIQUE

Les intégrales qui apparaissent dans le calcul des matrices élémentaires font intervenir des expressions compliquées. C'est pourquoi il est pratiquement exclu sauf dans quelques cas particuliers de les calculer analytiquement. Ces intégrales sont évaluées numériquement sous la forme

dtdn = Y f{t,r).Pr

( K 3 4 )

n

où E, r

sont les coordonnées d'un ensemble de points d'intégration sur l'élément réel (par exemple des points de Gauss)

Pr sont les coefficients d'intégration numérique (poids) associés à ces points.

O n notera que les tranformations effectuées et la technique d'intégration utilisée nécessite le calcul des expressions du type Nj (Ç r

), 9Nj(Ç r

)/3Ç, ... Ces expressions ne dépendent pas de la géométrie de l'élément réel et peuvent donc être calculées une fois pour toute pour chaque type d'élément.

7.5.4 - ORGANISATION DU CALCUL

Soit un élément isoparamétrique comportant n noeuds d'interpolation et q points d'intégration numérique. La première étape consiste à calculer la valeur des fonctions d'interpolation et de leurs dérivées aux q points d'intégration. Cette opération est réalisée par les sous-programmes :

G A U S S : stockage dans les tableaux V C P G et V K P G des poids et des coordonnées des points d'intégration sur l'élément de référence carré

G A U S S T : idem pour l'élément de référence triangulaire

14

F O N C - N I - 2 D : stockage dans le tableau V N I des valeurs

Le tableau V N I contient séquentiellement

<N

¡

>,<dN i

/dí,>,<dNJdr\> pour E, = E„ (1-34

1

)

< N > , < dN /di>, < dN / d n > pour E, = £, a v e c

< N . X U = < Ayç), A>

2

(£), - ,AMÇ) >

<ôiV ld£>(£) = O A T , (i,)\Bi,,...,dN (^)/8^>

(soit n x 3 x q valeurs)

Ce calcul n'est effectué qu'une fois.

La seconde étape consiste à calculer les quantités qui dépendent de la géométrie de l'élément réel. Ces opérations sont à effectuer pour chaque élément.

Dans le cas de la bibliothèque C V 2 D , les quantités à calculer sont : det (J(Ç r

)) (qui apparaît lors de la transformation des intégrales) et 3Ni(Ç r

)/3_

<

L e ca i eu i ¿es dérivées des fonctions d'interpolation sur l'élément réel fait lui-même intervenir le calcul du jacobien au point d'intégration. L'enchaînement des opérations est le suivant :

Pour chaque point d'intégration :

- calcul de la matrice jacobienne (2 x 2) à partir du tableau V N I et du tableau V C O R E (coordonnées des noeuds de l'élément) par la relation

(1.32)

- calcul du déterminant de [J] et de son inverse [J]~^. Ces opérations sont réalisées par le sous-programme J A C O B

- calcul (si nécessaire) des quantités 3Nj/3x et 3Nj/3y au point d'intégration. O n déduit de la relation (1.32)

ñN. ají.

771 _j r "«Ti

ñN i

\-

ÍJ]

[SNJ (1.34")

S y dn

Cette opération est réalisée par le sous-programme D N E D X à partir du tableau V N I et de l'inverse de la matrice jacobienne qui vient d'être calculée. Le résultat est stocké dans le tableau V N D C qui contient séquentiellement

15

< dN./dx > , < 8 N /dy>

i i

<dN./dx > , <dN./dy> pourí,= í

>

(1.34

m

)

(soit n x Z x q valeurs)

Le calcul d'une matrice élémentaire procède d'une façon analogue en bouclant sur les points d'intégration et en utilisant les tableau V N I et/ou V N I X ainsi que les tableaux contenant les constantes de milieu et de force de l'élément.

Cette organisation a été retenue car les tableaux à stocker (du type V N I X ) ne sont pas très importants. Il est ainsi possible d'utiliser le fait que certaines quantités telles que les jacobiens ont été calculées précédemment. C'est le cas lorsqu'on calcule dans la m ê m e boucle sur les éléments la matrice de rigidité puis la matrice de masse. C'est pourquoi les 10 premières variables du c o m m o n

/ L M T R A V / s o n t des indicateurs qui sont mis à zéro par le programme M E F I S T O lorsqu'on change d'élément. Il suffit alors de tester la valeur d'un indicateur pour savoir si le calcul des déterminants ou des jacobiens aux points d'intégration a déjà

été réalisé précédemment.

7.5.5 - CALCUL DE LA MATRICE DE RIGIDITE

L'équation effectivement traitée dans la bibliothèque C V 2 D suppose que la divergence de la vitesse est nulle sur le domaine soit V U = 0 . Cette hypothèse est vérifiée s'il n'y a pas de source de masse dans l'écoulement. Si les seules sources proviennent de puits ponctuels, cette hypothèse est encore vérifiée sur le domaine d'étude sauf aux noeuds qui représentent les puits (mais on sait que la vitesse n'est pas définie en ces points). Il s'avère donc que cette hypothèse n'est pas très restrictive en pratique.

Nous supposerons par ailleurs que le coefficient de diffusion moléculaire d o

et les dispersivités a^ et QVJ; sont constantes sur l'élément mais que la vitesse est variable. Suivant les conventions adoptées dans M E F I S T O et détaillées dans [1],

- d o

, otL et d j sont des constantes de milieu

- la vitesse est une propriété nodale dont il faut spécifier l'interpolation sur l'élément.

U n choix naturel consiste à adopter une interpolation de la vitesse identique à celle de l'inconnue variationnelle (c'est-à-dire la température ou la concentration). L a vitesse en un point de l'élément est alors calculée à partir de ses

16

valeurs aux noeuds de l'élément par :

U = U .N.

x

xj J

U =U .N

y yj j

Le c h a m p de vitesse ainsi défini est continu sur le domaine d'étude V . Il n'est cependant qu'une approximation du c h a m p de vitesse réel, ce qui est une des raisons pour laquelle l'équation traitée suppose initialement que div U = O .

Pour fixer les idées, prenons le cas du champ de vitesse induit dans un milieu homogène et isotrope par un ou plusieurs puits ponctuels. C e c h a m p est à divergence nulle (sauf sur les puits) mais son expression fait intervenir des fractions rationnelles qu'il n'est pas possible de représenter exactement avec des fonctions d'interpolation composées de polynômes. Il s'ensuit que la divergence calculée dans l'élément est susceptible d'être non nulle bien que le champ de vitesse original soit à divergence nulle.

Dans le cadre des problèmes susceptibles d'être traités pour l'aménagement, il est probable que le champ de vitesse ne sera pas connu de façon analytique. Plus précisément, il résultera d'un calcul numérique portant sur l'équation de la diffusivité auquel cas, c'est la valeur de la charge hydraulique (ou de la fonction de courant) qui sera connue aux noeuds du maillage (dans la mesure où le maillage hydraulique est identique au maillage utilisé pour résoudre l'équation de diffusion convection). L a vitesse en un point de l'élément s'obtient alors à partir des dérivées du potentiel hydraulique (et du tenseur de perméabilité dans l'élément). Si les fonctions d'interpolation utilisées pour le calcul hydraulique sont de classe C ° (ce qui est le cas des éléments de la bibliothèque D F 2 D cf [1]), le c h a m p de vitesse calculé est continu sur l'intérieur des éléments et discontinu à la traversée des frontières. O n notera cependant que le calcul de la matrice de rigidité est réalisé e n évaluant la fonction à intégrer aux points de Gauss et qu'il n'est donc nécessaire de connaître la valeur de la vitesse qu'en ces points. Par ailleurs, on peut montrer que les points optimaux pour calculer les vitesses à partir des dérivées des fonctions d'interpolation sont justement les points de Gauss et de fait, les algorithmes permettant de calculer un champ de vitesse continu sur le domaine s'appuient sur les valeurs calculées aux points d'intégration. Il s'ensuit qu'il est préférable de ne pas essayer de lisser les valeurs de la vitesse aux noeuds pour introduire un champ de vitesse continu mais plutôt d'utiliser directement les valeurs nodales préalablement calculées du potentiel hydraulique.

Le calcul de la matrice de rigidité élémentaire [r] = [k] + [z] (non

17 symétrique) s'effectue de la manière suivante :

- calcul préalable des det (J), < 8N¿/8x >, < 3Nj/3y > aux points de Gauss si nécessaire (cf. 1.5.4)

- W = [0]

- Boucle sur les points d'intégration

* Calcul de la vitesse (Ux, U y ) au point d'intégration (de coordonnées (Ç, ri)) ; 3options sont possibles :

a - on connaît la vitesse aux noeuds de l'élément, alors :

1

(1,34"")

.

I y

J

-

N

J i

^

)

b - o n connaît le potentiel hydraulique a u x noeuds et le tenseur d e perméabilité (constant) sur l'élément. O n suppose q u e l'interpolation d u potentiel h sur l'élément est identique à celle d e l'inconnue variationnelle soit*, h = y hN^

s

(ce qui est le cas si le calcul de h a été effectué avec la bibliothèque D F 2 D ) . O n a alors, d'après la loi de Darcy :

U -K y h. . dN ./dx + K y h . dN I dy

X XX ¿— J J IV — J J

J J

U -K y h . dN ./dx + K y h . dN ./dy

y

xy — j J yy "— J J

J J

où

KJQJ,

K

X y et Kyy sont les composantes distinctes du tenseur de perméabilité sur l'élément.

c - O n connaît la fonction de courant. O n suppose alors que ;p = ^_ ip . iV sur l'élément et que la vitesse est :

J

U - - dy/dy = - y ip .dN ./dy

x

• ~~

J J

en plan

U - d u i / d x = y w . d N ,/dx

y • —

J J j

18

en radial où r est le rayon (r = x ou y selon la

r

valeur de NAXIS dans le c o m m o n / P R O B / )

U y

= - -dylBx (cf[l]).

* Calcul du tenseur D au point d'intégration à l'aide de (l.Z).

* Modification de [r] pour i = 1... n et j = 1 ... n

r.. = r.. +

y y

dN./dx + {D ôN./dx + D dN./dy + f.U .N.)

I

XX 1 IV í

J

X l

où P^ est le poids du point d'intégration.

+ dN Jdy + (D dN./dx + D dN.id +f.U .N)

j

xy i yy ¡ >' y i

. dci ( u u , n)).

Pk

Le coefficient f qui intervient dans l'expression précédente permet d'annuler la vitesse sur un élément. Sa valeur normale (donnée dans les constantes de milieu MTT.T/SURF) est 1. Il permet de juxtaposer une région avec écoulement et une région sans écoulement tout en définissant la valeur de la vitesse aux noeuds. C e facteur doit être positionné à 0 dans une zone de diffusion pure.

1.5.6 - CALCUL DE LA MATRICE DE MASSE

Le coefficient w est supposé constant sur l'élément, n est possible de choisir entre deux formes de la matrice de masse :

- matrice de masse consistante ; c'est-à-dire celle qui est issue de l'application stricte de la méthode de Galerkine : r

m = w

N.N .dV = w Y N.d, ,n ).N (i, ,n. )det (J (i, , n )) . P

9=1

1

v

matrice de masse diagonalisée. Cette technique consiste à annuler les termes qui ne sont pas sur la diagonale mais il faut modifier les termes de la diagonale de manière à retrouver la masse totale de l'élément.

Cette simplification conduit à de meilleurs résultats pour certains types d'équations (équation de la diffusivité par exemple) mais n'est pas recommendée pour l'équation de diffusion-convection. Si on désire utiliser une matrice de masse diagonalisée, il faut positionner à 1 le paramètre N M D I A G du c o m m o n / P R O B / en appelant le bloc P R O B du

19 programme M E F I S T O . La matrice de masse calculée est alors : m.. = wN dV m = 0

il ! e i y

ce qui revient à affecter au coefficient m y la s o m m e des termes de la ligne i puisque

^ T A

7

. = 1 pour les fonctions d'interpolation utilisées.

j

1.5.7- CALCUL DES SOLLICITATIONS REPARTIES

Le terme de sollicitations réparties est supposé constant sur l'élément et son calcul est analogue à celui des coefficients diagonaux de la matrice de masse diagonalisée.

g q

% q q

9=1

1.5.8- CALCUL DES SOLLICITATIONS SUR LES FRONTIERES

Le flux h = DjjC,j.n^ est supposé variable sur la frontière. Les n^ étant les composantes de la normale extérieure au domaine, h est positif lorsque le flux diffusif est dirigé vers l'intérieur du domaine. Pour évaluer la contribution du flux sur la frontière, nous transformons d'abord l'intégrale qui est ensuite évaluée directement ou numériquement selon le cas.

- Transformation de l'intégrale

L'intégrale h.N..dS s'écrit, en fonction de l'abscisse curviligne sur l'arête de référence :

h.N ( Ç ) J

C

( U < U avec J

s

( U = Vidx/di,)

2

+

Nous considérerons ici le cas d'une arête linéaire ou quadratique.

2

il

g, i .

-1 0

arête de référence

20

La géométrie de l'arête est décrite par les fonctions d'interpolation de l'élément monodimensionnel linéaire ou quadratique selon le cas. Si l'arête est droite, le jacobien de la transformation géométrique est constant et vaut L / 2 , L étant la longueur de l'arête. Il est alors possible d'intégrer exactement la contribution du flux.

Si l'arête est courbe, le jacobien n'est plus constant sur l'arête et il faut intégrer numériquement (cf. sous-programme C V 2 D 7 Q de la bibliothèque C V 2 D dans lequel l'intégration est réalisée avec 3 points).

Dans le cas axisymétrique, on a dS = 21ir.dr.dz et l'intégrale à évaluer est

2 n h .N (i,) . r (Ç) J c

(£,) di

Nous donnons ci-dessous l'expression du vecteur f s

aux noeuds d'une arête droite.

- Arête linéaire

M

(le flux varie linéairement)

f

-

L

-

« 6

2 l

1 2

(en plan) f = 2n —

12 où rj est la valeur du rayon aux noeuds.

(cas axisymétrique)

- Arête quadratique droite h3 attention : h¿ est la valeur du flux au noeud 3 de l'arête

3 le flux varie de façon quadratique

h

== où les N ' k sont les fonctions d'interpolation de l'élément m o n o dimensionnel quadratique

* *

(en plan)

n

30

fs

L

~ 30

'

7

V

• r

2

(Tj + r

2> i

L.

r-

4

1

2

- 1

4

2

2

2

16

h z

L

h

2.

-

J

"V h v r

2

2r l

2 r

2

16 (r +

(en géométrie axisymétrique)

21

1.6 - UTILISATION DE LA BIBLIOTHEQUE CV2D

La structure de la bibliothèsque C V Z D est conforme à la description donnée dans [1], Elle comporte à l'heure actuelle 5 types d'éléments. Elle est composée des sous-programmes suivants :

E L E M L B (appel des différents types d'élément)

E L E M 0 1 - C V 2 D (exécution des fonctions élémentaires pour l'élément de type 1 : triangle linéaire 3 noeuds)

E L E M 0 2 - CVZD (exécution des fonctions élémentaires pour l'élément de type 2 : quadrangle linéaire 4 noeuds)

E L E M 0 3 - CVZD (exécution des fonctions élémentaires pour l'élément de type 3 : triangle quadratique 6 noeuds)

E L E M 0 4 - CVZD (exécution des fonctions élémentaires pour l'élément de type 4 : quadrangle quadratique incomplet 8 noeuds)

E L E M 0 5 - CVZD (exécution des fonctions élémentaires pour l'élément de type 5 : quadrangle quadratique Lagrange 9 noeuds)

C V 2 D 3 (calcul de la matrice de rigidité pour un des éléments de la bibliothèque)

C V 2 D 5 (calcul de la matrice de masse)

C V 2 D 6 (calcul du résidu [r] {C})

C V 2 D 7 (calcul du second m e m b r e Q . N . dV

r ¡v

l

C V 2 D 7 L (calcul de h.N .dS

s u r u n e a r

®te linéaire)

•'S î

CV2D7Q (calcul de sur une arête quadratique)

C V 2 D 8 (calcul et impression les gradients aux points d'intégration)

C V 2 D 9 (calcule les vitesses aux points de Gauss et les stocks dans le tableau des propriétés élémentaires pour utilisation future par C V 2 D 3 et

CV2D6.

C V 2 D K X (calcul du déterminant du jacobien et des dérivées des fonctions d'interpolation aux points d'intégration).

Ces sous-programmes se trouvent dans la bibliothèque L M C V 2 D . O L B .

Les sous-programmes généraux de MEFISTO se trouvent dans la bibliothèque M E F . O L B .

Le programme principal (qui réalise l'appel des différents blocs fonctionnels de MEFISTO) est P P M E F . O B J .

22

La création du module exécutable M E F . C V 2 D est réalisée par l'instruction

L D J K / E X E = M E F C V 2 D P P M E F , L M C V Z D / I N C = E L E M L B / L , M E F / L , L I B C O M / L .

On trouvera ci-dessous l'entête du sous-programme E L E M L B de la bibliothèque C V 2 D qui précise la façon de définir les paramètres de l'équation à l'aide des blocs M U I , F O R C , P R N D .

23

2 - NOTIONS SUR LES SCHEMAS UPWIND, LA STABILITE ET LA DISPERSION

DES SCHEMAS NUMERIQUES

Les solutions numériques de l'équation de diffusion-convection présentent souvent des oscillations lorsqu'elles sont calculées par la méthode des éléments finis avec la formulation de Galerkine (que nous appellerons M E F G dans la suite). Des problèmes identiques apparaissent avec la méthode des différences finies centrées

( M D F C ) . Les causes de ces oscillations sont diverses et nous nous proposons d'indiquer brièvement les moyens d'y remédier.

2.1 - LES SCHEMAS "UPWIND"

2.1.1- CAS MONODIMENSIONNEL

Pour fixer les idées, nous étudierons rapidement la discrétisation de l'équation de diffusion-convection stationnaire monodimensionnelle à coefficients constants, soit :

dT d

2

T

U = D 0 < x < L (U>0)

dx

dx~ (2.1) avec les conditions aux limites suivantes

Ti0) = T

° (2.2)

T(x = L) = T

L

et T

L

*T

o ( 2 > 3 )

La solution exacte de l'équation précédente est :

p

x

-

T (x)-T i e

*

L

T* = °- = — (2.4)

T - T

p

L

» 1 - e

e

où P e

= U . L / D est le nombre de Péclet qui mesure l'importance du transport convectif par rapport au transport diffusif.

Tant que P e

est faible, la solution varie "doucement" sur le domaine mais lorsque P e

croit, la solution devient du type couche limite, c'est-à-dire T(x) = T o

sauf au voisinage de l'extrémité x = L où elle s'adapte brutalement à la condition d'exhaure

T (L) = T L (fig. 2.1)*. L'épaisseur de cette couche limite peut-être estimée en égalant les termes convectif s et diffusif s de l'équation sur une faible distance 6, soit :

U.àT DAT

* En fin de texte

24 ce qui conduit à 5/L = 1 / P e

qui indique que l'épaisseur ô de la couche limite est approximativement inversement proportionnelle au nombre de Péclet.

La difficulté de ce problème réside dans le fait que toute la physique de la solution est contenue à l'intérieur de cette mince couche limite lorsque P e

>> 1.

• La discrétisation de l'équation (2.1) par la M D F C ou la M E F G linéaire (i.e.

avec des éléments linéaires) sur une grille uniforme de pas h donne au noeud i : soit:

2 h

P

— (T. - T. , ) = T , - 2 T + T. , (2.5)

2 i + i i - i i - i i i + i

U.h (h

P

est le nombre de Péclet de maille et N le nombre d'éléments de la grille.

La solution exacte de l'équation discrétisée (2.5) avec les conditions aux limites associées (2.2) et (2.3) est :

I

T i

- T

- T

0

0

1 -1

(-

+ Pn y

1

! (

\

"I

- 9,

+_Pn_yv

-Pn )

La forme de cette équation montre clairement que pour P n

> 2, des oscillations qui n'existent pas dans la solution réelle (2.4) vont apparaître. L'aspect de ces oscillations dépend de la parité de N et leur amplitude de la valeur P n

.

Il existe de nombreux schémas de discrétisation appelés "upwind" permettant de supprimer les oscillations lorsque P n

est > 2. C e s schémas partent du principe que c'est la valeur de la variable en amont qui va conditionner la valeur au noeud de calcul (et ceci d'autant plus que la convection sera importante) et consistent

à décentrer vers l'amont le terme convectif. Ainsi, le schéma "upwind" pur s'écrit

(pour u > 0) :

Pn(T.-T -) = T. - 2 T + T (2.8)

I I - 1

l

- 1 1 i + 1

Les meilleurs schémas upwind sont basés sur des équations discrétisées dont la solution satisfait toujours l'équation (2.4) aux noeuds. Ils sont obtenus, dans le

25

cas d'une discrétisation en éléments finis linéaires, en utilisant des fonctions de pondération quadratiques asymétriques ou en déplaçant la position des points d'intégration numériques sur l'élément de référence. Quelle que soit la technique utilisée, ces schémas optima se présentent sous la forme (identique pour la M D F C et la M E F G linéaire).

Pn \ tanh — )(T_,-T

:

, ) = T. , - 2 T . . +

2 / 1 + 1 l - T 1 - 1 ! 1 + 1

(2.9)

Cette équation est identique à (2.5) si ce n'est que le nombre de Péclet de maille effectif P ^ est :

K=2.tanh{PrU2)

( 2 > 1 Q )

Si P n

est remplacé par P n

dans l'équation (2.5), la solution obtenue aux noeuds sera identique à la solution exacte (2.4). Le schéma (2.9) exact (dans ce cas simple) s'applique aussi au cas d'une grille à pas variable sous la forme :

T - T =coth(Pn /2) (T , , - T ) + cothiPn 12) (T .-T.) , „ .

i + i i i + i i + i i i i - i i ( 2 . 1 1 )

où P n i = U.hj/D est le nombre de Péclet de maille sur l'élément i.

L'examen de la relation (2.10) montre que P n

= P n

lorsque P n

<< 1 et P n

=

2 lorsque P n

>> 1. Par ailleurs, la diffusivité originelle D de l'équation (2.1) est remplacée par la diffusivité effective

_ Pn Pn

D = D . — coth —

2 2 (2.12) qui donne D -> D lorsque P n

<< 1 et D -> 1/2 U . h lorsque P n

>> 1. O n voit que lorsque

P n

est grand (P n

> 5 environ) la diffusivité effective devient indépendante du nombre de Péclet.

Le gros inconvénient de ces schémas upwind est d'introduire une diffusion artificielle (diffusion numérique) qui devient largement supérieure à la diffusion physique lorsque l'écoulement est dominé par la convection.

Pour conclure cet exemple, les oscillations introduites par un schéma centré lorsque P n

> 2 peuvent être interprétées c o m m e un avertissement. La solution devient difficile dans la zone près de l'exhaure parce qu'il s'y développe une couche limite dont l'épaisseur est faible par rapport au pas de la grille utilisée. Cet avertissement n'est pas donné par un schéma upwind et l'utilisateur non averti peut

26

être a m e n é à croire à la validité de la solution pour n'importe quelle valeur du nombre de Pédet. Le remède aux oscillations consiste à suivre correctement la couche limite lorsque P e

est grand. Ceci peut être facilement réalisé en affinant (au moins grossièrement pour limiter l'amplitude des oscillations) le maillage au voisinage de l'exhaure. D e fait, dans le cas considéré, toute la solution "réside" dans la couche limite et elle est pratiquement triviale à l'extérieur de la couche limite. Il est alors logique de positionner des noeuds de calcul dans la région contenant la couche limite.

2.1.2 - EXTENSION A 2 OU 3 DIMENSIONS

II est clair que la comparaison des schémas pour résoudre un problème monodimensionnel est quelque peu stérile si ces schémas et leurs résultats ne sont pas généralisable à 2 ou 3 dimensions. Dans le cas de la diffusion-convection stationnaire, il apparaît que des généralisations sont souvent possibles bien que la mise au point de schémas upwind multidimensionnels vraiment efficaces soit plus difficile.

La difficulté rencontrée dans l'exemple monodimensionnel (condition de

Dirichlet sur un exhaure associée à un nombre de Pédet élevé) apparaît aussi dans les cas Z D ou 3 D . Heureusement, il est rarement nécessaire en pratique d'imposer ce type de condition sur une limite aval. En effet, l'exhaure est souvent une limite fictive qui est introduite parce qu'il n'est pas possible de mailler le domaine d'étude qui est infini. L a véritable condition à la limite est alors inconnue. L'expérience montre que la "meilleure" condition à imposer sur ce type de frontière est celle qui est naturellement introduite dans la formulation variationnelle (i.e. flux diffusif normal nul). Cette condition de flux nul, dont l'efficacité sera montrée plus tard sur un exemple, ne génère pas de couche limite au voisinage de l'exhaure et apparaît plus réaliste qu'une condition de Dirichlet. Dans les quelques cas où il serait plus réaliste d'imposer la valeur de la concentration ou de la température sur l'exhaure, l'utilisateur doit savoir qu'une couche limite importante risque d'apparaître et qu'il devra donc discrétiser cette zone en conséquence.

Les schémas upwind bidimensionnels introduisent souvent une diffusion numérique très importante (qui rend de fait les résultats indépendants du nombre de

Pédet dès que la convection domine). L'effet est généralement encore plus catastrophique lorsque la vitesse n'est pas parallèle aux lignes du maillage avec l'apparition d'une diffusion numérique transversale à l'écoulement. O n notera cependant qu'il existe certains schémas (utilisant notamment la méthode des caractéristiques ou une approximation spatiale plus précise du terme convectif) qui

27

sont très efficaces dans certaines situations.

2.2 - NOTIONS D E STABILITE, DISPERSION ET DISSIPATION DES S C H E M A S

NUMERIQUES

Les schémas upwind brièvement présentés dans le paragraphe précédent ont généralement été développés pour résoudre l'équation de diffusion-convection stationnaire avec des nombres de Péclet élevés. L'intégration temporelle des

équations discrétisées dans l'espace introduit de nouvelles approximations. D e nombreux travaux sont actuellement consacrés à l'étude de la stabilité, de la dispersion et de la dissipation des schémas numériques. Nous nous proposons d'examiner ici deux exemples très simples pour introduire ces caractéristiques des schémas numériques.

2.2.1 • STABILITE (EQUATION DE LA DIFFUSION)

L'analyse des caractéristiques d'un schéma considère généralement que les discrétisations spatiales et temporelles sont effectuées indépendamment l'une de l'autre. En règle générale, on réalise d'abord la discrétisation spatiale ce qui conduit à un système d'équations différentielles ordinaires qu'il faut intégrer. Bien que le but soit de résoudre le système d'équations aux dérivées partielles originel, la majorité des analyses s'intéresse aux propriétés de l'opérateur d'intégration temporelle, certaines caractéristiques du système discrétisé dans l'espace étant connues et donc à la solution de types particuliers de systèmes d'équations différentielles. H importe cependant de garder à l'esprit que toutes les caractéristiques de la discrétisation

(spatiale et temporelle) doivent normalement être considérées simultanément.

Pour illustrer le concept de stabilité, nous nous basons sur le système d'équations différentielles ordinaires suivant :

KT=F (2.13) dans lequel la matrice M (matrice de masse) est symétrique définie positive et la matrice K (matrice de rigidité) est semi-définie positive. L a discrétisation spatiale de l'équation de la diffusion conduit à un système d'équations présentant ces caractéristiques.

28

Le schéma d'intégration temporelle utilisé dans les blocs T R L C et T R L V du programme M E F I S T O est le suivant :

(M+BAtK) AT= -KT l

+ F(t i

+ 9 Ai)

{ Z

'

1 4 ) où At est le pas de temps

6 est la pondération implicite-explicite du schéma d'Euler.

Quelques valeurs particulières de 6 sont données ci-dessous :

9 = 0 schéma explicite

6 = 1 schéma implicite pur

9 = 1/2 schéma de Crank-Nicholson

9 = 2/3 schéma de Galerkine (qu'on obtient en discrétisant le système (2.13) avec des éléments finis linéaires)

L'influence de 8 peut-être étudiée en transformant le système (2.13) en un système découplé grâce à la transformation

T = X V

(2.15) où la matrice de transformation X est constituée par les n vecteurs propres (n : nombre d'équations du système (2.13)) Xj définis par :

(K-\.M)X. = 0 ¿ = l , . . . , n (

2

-

16

)

La matrice X est orthonormée et vérifie

X

(

M X =8

i j ij

X '. K X j

- \ x

5.. (pas de sommation) (2.17)

Dans la base des vecteurs propres, le système (2.13) s'explicite sous la forme des n équations découplées :

V k +

h

V

k = *k <

2

-

18

>

On peut alors examiner l'influence de 6 en étudiant l'équation simple

29

(2.19) ce qui permet de comprendre le comportement des différents modes qui forment la solution du problème initial.

L a solution analytique de (2.19) est :

L'équation (2.19) discrétisée par le s c h é m a d'Euler donne

( i + e A A ¿ ) v

I + 1

= ( i - ( i _ e ) A A O v ' soit V

l + l

= A . V l

(2.21) où i _ ( i _ e ) A At (2.22)

A = i + e

A

At

est appelé facteur d'amplification.

D a n s le cas où le s c h é m a d'intégration temporelle fait intervenir la valeur de la variable à plus d e deux instants, A est une matrice (matrice d'amplification).

L'étude de la stabilité, de la dispersion... du s c h é m a est alors basée sur les valeurs propres de la matrice d'amplification.

L a solution exacte (2.20) indique que

\ V l + l

\ < \ V l

\ si A > 0

y i + i _ y '

SL

x = o (2.23)

(compte tenu des caractéristiques du système (2.13), les modes X sont réels et > 0).

Le schéma sera donc stable si IA | < 1 (la deuxième condition est satisfaite si X = 0) soit i _ d _ e i

A A<

- 1 < < 1

1 + 6 X 4 1

UM)

L'inégalité de droite est toujours satisfaite pour toutes les valeurs possibles des paramètres (X At > 0 et 0 s 6 < 1) et l'inégalité d e gauche est vérifiée si 6 > 0 . 5 . Si 6

< 0.5, l'inégalité de gauche impose

30

Xàt < 2/(1 - 2 6) (8 < 1/2)

(2.25) ce qui, pour X donné, fixe une valeur maximale du pas de temps.

Le schéma d'Euler est dit conditionnellement stable pour 6 < 1/2 et inconditionnellement stable pour 6 > 1/2.

R e m a r q u e : II est impératif de vérifier la condition | A l < 1. En effet, la solution de

(2.21) peut s'écrire V * = A * V o

et si I A | est > 1, toute erreur (résultant ne serait ce que de la précision limitée d'un ordinateur) va croitre indéfiniment.

Dans le cas conditionnellement stable, la condition de stabilité At < 2/((l-

29)X) doit être vérifiée pour tous les m o d e s (i.e., tous les Xjç, k = 1, ..., n) du système

(2.13). L e plus grand m o d e ( X m a x

) impose la restriction la plus sévère sur le pas de temps maximal. E n fait, dans le cas de l'équation de la diffusion, on peut montrer que

X m a x

= 0(h~2) (i.e. est de l'ordre de h~¿) où h mesure le pas de la discrétisation spatiale. L e pas de temps m a x i m u m doit donc vérifier At < constante, h ^ . Si le système comporte beaucoup d'équations, cette restriction est sévère et c'est pourquoi les algorithmes inconditionnellement stables sont généralement utilisés.

R e m a r q u e : Si le système (2.13) est issu d'une discrétisation par élément finis, alors

X m a x

est inférieur à la valeur propre maximale des éléments pris individuellement.

Le comportement de A en fonction de XAt est représenté sur la figure 2.2

pour différentes valeurs de 6 . L a valeur de XAt pour laquelle A = 0 est appelée la limite d'oscillation parce que pour des valeurs supérieures, le signe de A

1

change à chaque pas. Pour les algorithmes inconditionnellement stables (6 ¿ 1/2) la valeur asymptotique du facteur d'amplification, A « , vérifie I A I » s 1. Si 8 > 1/2, alors quel que soit XAt, IA | < 1 et les composantes modales élevées (qui sont généralement erronées et donc indésirables dans les schémas numériques) décroissent. Cependant si

6 = 1/2 (ou est très près de 1/2) et XAt >> 1, alors A = -1 et les composantes modales

élevées se comportent c o m m e (-1)

1

. C e s oscillations dans le temps se manifestent souvent dans les calculs. O n peut dans ce cas filtrer ces composantes élevées en utilisant la moyenne de deux instants ( V i + 1

+ V^/Z puisque (-1) + (-l) i+1

= 0.

R e m a r q u e : Si a = 1 et At -> °°, on obtient la solution stationnaire (i.e. celle pour laquelle V = 0).

31

Le schéma d'Euler est consistant, c'est-à-dire que l'erreur de troncature locale (qu'on obtient en remplaçant V* et V*

+ l par leur valeur exacte dans l'équation discrétisée) vérifie I e(t) I < C . At^

+

^ pour tout t. C est une constante indépendante de At et k > 0 est l'ordre de précision (ou taux de convergence). O n montre par ailleurs facilement que le schéma d'Euler vérifie k = 1 pour tout 8 e [0,1] sauf pour

6 = 1/2 auquel cas k = 2. Cette propriété associée à la stabilité est une condition nécessaire et suffisante qui assure la convergence du schéma (i.e. l'erreur tend vers zéro lorsque At tend vers zéro) (théorème de Lax).

2.2.2 - DISSIPATION, DISPERSION

Cette section se propose d'introduire les concepts de base associés à la dispersion et à la dissipation des schémas numériques sur un exemple simple. E n fait, la dispersion est liée à un type de solution plutôt qu'à un type d'équation. Plus précisément, un système dispersif est un système qui admet des solutions du type

T=Ae ilk

*wn

(2.26)

La fréquence U) est une fonction réelle du nombre d'onde k, A est l'amplitude de l'onde et : = r M

f

£a vitesse de phase (qui est différente de la vitesse d'onde en général) est c = 7 (2.27) et les ondes sont dites dispersives si c dépend de k.

La solution résulte en général de la superposition de plusieurs modes du type (2.26), ou, dans le cas le plus général, d'une intégrale de Fourier. Si la vitesse de phase c n'est pas la m ê m e pour tous les nombres d'onde k, les modes seront propagés à des vitesses différentes, d'où le phénomène de dispersion. L a quantité kx-u)t est la phase, 2 H / k la longueur d'onde et 21i/u) la période.

Pour un système non dispersif, la vitesse de phase est une constante C o

et w = c o k (la solution (2.26) est alors celle de l'équation monodimensionnelle hyperbolique). O n peut alors adopter la définition classique : la fonction (2.26) est dispersive si d^u/dk^ ^ 0.

Dans la plupart des cas, l'équation à traiter est discrétisée par rapport aux variables spatiales et temporelles. Chacune des discrétisations introduit de la dispersion et dans l'étude de la qualité de la solution numérique, c'est la dispersion

32 totale qui doit être considérée. Pour fixer les idées, nous considérons l'équation de diffusion-convection monodimensionnelle homogène (sans terme source) où U et D sont constants.

dT dT â

2

T .

— + U — - D = 0 Z.22,)

dt dx

Ô X

2

C e système est non dispersif. Pour une solution initiale de la forme :

U

'

0) = Te

(

2

.

29) la fréquence associée est u) = k U et le paramètre d'amortissement est £ = D k ¿ .

Nous considérons deux discrétisations spatiales de l'équation sur une grille

à pas constant h avec des éléments finis linéaires et la formulation de Galerkine est l'équation au noeud j qu'on obtient en appliquant directement la méthode. Si on diagonalise la matrice de masse (en affectant à un terme diagonal la s o m m e des termes de son rang) ^-, ^— f f f

m . = > m . = > \ N. N . dV = \ N Z N dV = N dV l'équation précédente s' l

' ' y

!

-> ^ I '

J

\ > ¿ - J J » l'équation précédente s'écrit :

T.+ — ( T ,-T. , ) - — ( T . -IT +T ,) = 0

J

2h J+l J ~

l h

2 J - i J j + l

Cette équation est identique à celle obtenue par les différences finies centrées. Les deux équations précédentes peuvent se mettre sous la forme :

( 1 + r L )

^

+

^

( r

,

+ 1

r

, -

1 )

- ^

L r

,

= 0

(2.30) où L est l'opérateur défini par L T ; = T;_i - ZTj + Tj+i, r = 0 pour la cas M D (matrice masse diagonale) et r = 1/6 pour le cas M C (matrice masse consistente, c'est-à-dire issue de l'application consistente de la formulation de Galerkine).

O n suppose que la solution de (2.30) peut se mettre sous la forme

W *

( t )

(2.31)

33

En injectant (2.31) dans (Z.30) on obtient r L ) S . . $ ( i ) + - ( S .

+ 1

- S j

_

1

) $ - soit :

— (S.

-S .)- —L S

2 / 1 ;

+ 1

J-

1

h

2 J

(1 + rL) S puisque le terme de gauche ne dépend pas de t, ce qui donne les deux expressions suivantes :

U

ß(l + rL) S + —

(T.

J 2h J

D

-T . ) - —L S =0

J

1

2

J

= 0

(2.32)

(2.33)

O n a :

Injectons la condition initiale (2.29) dans (2.32) pour calculer

A -

ikx

. >•»

L S =T Le

J

= T

,k( x . - h) i kx . i k{ x

+

h

)

J

- 2 e

J

+ e

J

= Te

=

Te

2 cos k h - 2

S

x l

- S

ikx.

J

ik h - ikh

e — e s* ikx -

= 2 ¿Te

J

sin k h

O n en déduit que ß est donné par l'expression

ß = -

lUsinkhlh + 2D (1 - coskh)/h

2

1 + 2r (coskh - 1 )

ß = -

iwsinkh/kh + 21, (1 - coskh)l(kh'

1 + 2 r (cosfc/i - 1 ) soit

(2.34) où u = ku et Ç = D k ^ sont la fréquence et le paramètre d'amortissement exacts vus précédemment.

La solution exacte de l'équation (2.33) est

(2.35)

34 et par analogie avec la solution de l'équation initiale (2.28) où

(2.36) nous récrivons l'équation (2.35) sous la forme :

< î > U ) = e lW

' (2.37)

II résulte de l'équation (2.34) que le rapport de la fréquence semidiscrétisée UJ à la fréquence exacte w est :

w sin khi kh

w ~ 1 +2r(cos kh- 1 ) (2.38) et que le rapport du paramètre d'amortissement semi-discrétisé ç au paramètre d'amortissement exact ç est

l 2 ( 1 - cos k h ) / (kh )

2

l~ 1 + 2 r(cos kh - 1 ) (2.39)

Le rapport des fréquences (2.38) et le rapport d'amortissement (2.39) sont fonctions du n o m b r e d'onde adimensionnel kh. L e résultat important qui apparaît sur la figure 2.3 est que les deux discrétisations ( M D et M C ) réduisent les fréquences.

Le s c h é m a M D exhibe une réponse en fréquence particulièrement mauvaise et on doit donc s'attendre à obtenir une mauvaise solution numérique. Pour ce qui concerne le paramètre d'amortissement, le schéma M C l'augmente alors que le s c h é m a M D le réduit (fig. 2.4).

L'application d'un opérateur de discrétisation temporelle conduit à discrétiser les équations semi-discrétisées (2.32 et 2.33). Nous appliquons le s c h é m a d'Euler à l'équation (2.33) où ß = - (i u) + ç) ce qui donne

4> ( t + dt ) 1 - ( 1 - 6 ) Ai ( i w + Û

= = (2 40)

4>(i) 1 +Bàt (i w + l)

Par analogie avec la solution exacte de (2.33) qui est donnée par (2.37), nous réécrivons 4>(t + dt) = e~^ w +

^<J)(t) ce qui donne, c o m p t e tenu de (2.40)

e

1 + 9 Ai (i w + Q dont on déduit que

tg w Ht - b I a

(2.41)

35

et e - c, Ai = Va

2

+ b

2

1 c

(2.42) avec a = [l -(l -Q)Càt][l + QlAt]-e(l -d) (Z M )

2

(2.43)

b = w A T (2.44) c = ( i + e < Ai)

2

+ e

2

( w A o

2 ( 2 # 4 5 )

Le rapport u / œ dans le cas de la convection pure (ç = 0) est représenté sur la figure 2.5 en fonction de coAt pour différentes valeur de 6. O n remarque qu'un schéma implicite avec 6 > 0.5 réduit toujours la fréquence semi-discrétisée w .

C o m m e la discrétisation spatiale réduit aussi les fréquences (fig. 2.3), il faut utiliser le schéma avec la matrice masse consistente (MC) pour limiter la dispersion lorsqu'on utilise le schéma d'Euler avec 9 > 0.5.

Remarque : l'expression (2.42) montre que l'intégration explicite (9 = O ) dans le cas de la convection pure (D = 0) conduit à un schéma inconditionnellement instable. O n trouve en effet que dans ce cas : alors que le schéma ne peut être stable que si D > 0.

Une mesure de la dissipation numérique introduite par la discrétisation temporelle est donnée par le rapport Ç/UJ. O n remarquera sur la figure 2.6 que le schéma de Crank-Nicholson n'amortit pas les fréquences semi-discrétisées dans le cas de la convection pure.

Les résultats concernant la dispersion introduite par la discrétisation spatiale (mesurée par les variations du rapport û / w en fonction de kh) s'étendent qualitativement au cas bidimensionnel. L a dispersion introduite par un schéma dans lequel la matrice masse est diagonalisée est toujours beaucoup plus sévère que dans le cas où la matrice de masse consistante est utilisée. Cet effet est particulièrement désastreux dans le cas des quadrilatères bilinéaires (4 noeuds) et biquadratiques

incomplets (8 noeuds) (cf. [2]).

2.3 - OSCILLATIONS A V E C L ' E Q U A T I O N D E LA DIFFUSION

Nous reprenons ici un exemple donné dans [3] qui montre que la méthode de Galerkine "consistante" peut produire des oscillations dans le cas de l'équation de la diffusion. L'exemple présenté montrera que ces oscillations peuvent être interprétés c o m m e un avertissement.

36

Considérons le domaine (0 < x < 1) discrétisé sur une grille à pas constant h formée d'élément linéaires et le problème :

a

2

T

(2.46)

dt

T = 0 àt = O

T = 1 à x = O

T = 0 à x = 1

c'est-à-dire que les conditions aux limites se traduisent par une variation brutale de la température à x = 0 à l'instant t = 0

+

. Pour l'équation de la diffusion et un élément linéaire, les matrices de masse et de rigidité élémentaires sont :

m e

}= ~

\m e

]= -

2 1

1 2

1 0

0 1

dans le cas M C (matrice masse consistante) (2.47) dans le cas M D (matrice masse diagonale) (2.48)

1 - 1

- 1 - 1

où D est le coefficient de diffusion

(D = 1 dans le problème traité)

(2.49)

Si la grille de calcul comporte N + 1 éléments (de taille h = 1/N + 1 ) , l'assemblage des matrices élémentaires conduit, après élimination des équations correspondant aux degrés de libertés bloqués en x = 0 et x = 1, au système de N

équations différentielles

MT= -KT + F

/ •

(2.50)

dans le cas d'un schéma M C (où T } est la température du 1er d.i. libre, i.e. à x = 1/h.

Dans le cas du schéma M D , on a M = hi ou I est la matrice identité (la discrétisation M D est dans ce cas strictmeent identique au schéma différences finies pusique les seules conditions aux limites sont du type Dirichlet).

O n peut montrer (cf. [3] pour les références) que l'inverse de la matrice masse consistante possède les propriétés suivantes :

37

- ses éléments diagonaux sont positifs

- les éléments non diagonaux décroissent en valeur absolue (dans le rapport ^.268) en fonction de leur distance à la diagonale et surtout, leurs signes oscillent.

L a valeur de T (t = 0) est donnée par M ~ " . C o m m e seul le premier

élément de f est non nul, T(t = 0) est égal à 1/h fois la première colonne de M ~ l , soit

"I; (t = 0) = M j ^ / h . C o m p t e tenu des propriétés d e M ~ l rappelés ci-dessus, TÎ (0) est > 0 si j est impair et < 0 si j est pair, l'amplitude de Tj (0) décroissant avec j croissant. Il s'ensuit que la température à tous les noeuds pairs "démarre" dans la mauvaise direction puisque la solution exacte de (2.46) ne fait pas apparaître de valeurs négatives de T (x, t).

D a n s le cas du s c h é m a M D (matrice masse diagonalisée, on a Tj (O) = et ij (0) = 0 , j = 2 , . . . , N c'est-à-dire qu'aucune des températures nodales ne démarre dans la mauvaise direction. Il semble donc à première vue que la m é t h o d e de

Galerkine consistante prête à caution.

C o m p a r o n s la solution exacte de (2.46) aux solutions numériques obtenues en résolvant de façon exacte le système (2.50) dans les cas M D et M C .

La solution exacte de (2.46) est

r>

œ

• 2 9

T ix t) = ( 1 - ) - - y

Sm nnX

- " n~ f

La solution exacte du système (2.50) est

. N

T (t) = ( 1 — ) — > sin j na e ¡-y ro\

J N+l N+ 1 ^ - , 1 - cosna (2.5ZJ

n = 1

où a = n / ( N + 1 ) et les valeurs propres de K ~ * M sont

MC 9 1 - cos n a

X n

=6(N+l) pour le schéma M C (2.53)

£ r cos n Q

X

MD

= 2 ( N + 1 )

2

1 - cos n a

p o u r

i e s c h é m a M D (2.54)

R e m a r q u e : o n peut retrouver les valeurs propres d u système semi-discrétisé à partir du rapport ç/Ç calculé p r é c é d e m m e n t (cf. 2.2.2). O n a en effet trouvé :

l 2(1 - coskh)/(kh)

2

l ~ 1 + 2r(coskh - 1) qui ne dépend pas de la vitesse de convection U .

38

L'expression (2.51) montre que les nombres d'onde de la solution exacte sont niT (n = 1, ..., <») soit k = n1i dans l'expression (Z.55). L e paramètre d'amortissement exact est ç = D k

2

(soit n^ 112 ¿ a n s

\

e c a s

présent). E n tenant compte du fait que kh = n1l/N + 1 = n o et h = 1/N + 1, on trouve

1=2 (2V+ 1)

2

(1 - cosnaVa + 2 r(cosna - 1)) qui est identique à (2.53) et (2.54) avec r = 1/6 et r = 0 respectivement.

Il s'avère donc que le schéma M C amplifie le paramètre d'amortissement exact (n~ f ¿) et que le schéma M D le réduit (dans un rapport sensiblement identique).

D e fait, la valeur moyenne de X

n

et X

n

est une bien meilleure approximation de n

2

Les résultats pour N = 4 (5 éléments de dimension 0.2) sont présentés sur les figures 2.7 et 2.8 qui sont extraites de [3]. O n observe ce qui a été prévu précédemment, à savoir que dans le schéma M C , les températures aux noeuds 2 et 4 sont négatives dans un premier temps. La figure 2.7 montre qu'au temps t = 0.004, la solution M C présente encore des oscillations alors que la solution M D paraît

"raisonnable". A partir de t - 0.02 (d'après la figure 2.7), toutes les températures de la solution M C sont positives et le restent (bien que la solution M D paraisse un peu plus précise).

Nous suivrons l'opinion de Gresho et al. [3] en disant que les oscillations du schéma M C pour les temps faibles indiquent que pour ce problème et ce maillage, il existe un laps de temps minimum (= 0.02 dans ce cas) pendant lequel la solution discrétisée est trop imprécise pour être utile. O n observera en effet sur la figure 2.8

que la solution M D est très mauvaise à t = 0.004 bien qu'elle ne présente pas d'oscillations et puisse donc sembler raisonnable.

2.4 - LE PRINCIPE DU MAXIMUM

2.4.1 - DEFINITION

Les valeurs négatives des températures nodales obtenues par le schéma

M C indiquent que ce schéma ne vérifie pas le principe du m a x i m u m que vérifie la solution exacte. Supposons par exemple que la variable soit la concentration d'une substance dans le domaine d'étude et qu'il n'y ait pas de source (ou de puits) de cette substance. Le principe du m a x i m u m énonce alors que la concentration m a x i m u m apparaît sur la frontière ou à l'instant initial et que la densité ne devient jamais

39 négative. Si le schéma de discrétisation ne vérifie pas le principe du m a x i m u m discret, alors la solution calculée peut prendre des valeurs négatives et osciller dans le temps et dans l'espace.

Les solutions de l'équation de la diffusion-convection

dT

— + UVT-DAT= Q surV

ôt dT/dn = 0 sur S

Q

T = g sur S c

T(x,0)-T surV

o

(2.56) vérifient aussi le principe du m a x i m u m qui s'écrit :

max max-T.max-g

V o

+ t max

0 , max Q

min \ min— T ,min„g

V o

t miniO, min Q ) s T ( x , 0 )

avec E = Sq x [0, T].

(2.57)

2.4.2 - CAS DE L'EQUATION DE LA DIFFUSION

O n peut montrer facilement que l'élément linéaire monodimensionnel vérifie le principe du m a x i m u m discret pour l'équation de diffusion stationnaire. Il faut pour cela utiliser le théorème suivant : ([4]) :

Si la matrice réelle inversible d'ordre n A = (a^;) est à diagonale dominante avec ay < 0 si i f- j et &a > 0 si i e [1, n] alors A ~ * > 0 (une matrice n x n est dite à diagonale dominante si

a s 0 Vi ).

La matrice de rigidité élémentaire de l'élément linéaire monodimensionnel pour l'équation de la diffusion s'écrit :

D

1 - 1

- 1 1

(2.58) ou h est la longueur de l'élément et D le coefficient de diffusion.

Il est clair que la matrice de rigidité globale (k) résultant de l'assemblage des matrices élémentaires vérifie le théorème précédent et donc que T = K " ^ . F > 0 si

F > 0.

40

Les éléments d'ordre supérieur ne vérifient pas le principe du m a x i m u m discret. E n effet, les matrices de rigidité élémentaires de ces éléments sont :

*• =-2-

1 - 8

7 - 8

- 8 16 pour l'élément quadratique

/ 3 2

(2.59)

et

120/1

AAA - 39

y m

pour l'élément cubique

- 5fi7 162

162 - 567

1296 - 891

1296

(2.60)

La présence de termes non diagonaux > 0 dans ces matrices élémentaires va entraîner la présence de termes non diagonaux > 0 dans la matrice globale.

/TI .

Dans le cas transitoire, le principe du m a x i m u m discret sera vérifié si la matrice M + a A t K vérifie le théorème p r é c é d e m m e n t énoncé. Les matrices de masse

élémentaires consistantes et diagonalisées des éléments linéaires, quadratiques et cubiques sont :

h

6

1 *?

1 O

O 1

élément linéaire

(2.61)

'C

m.

h_

30

4

- 1 - 2

4 2

16

h

6

élément quadratique

(2.62)

h

1680

128 19

128

sym

99

- 36

648

- 36

99

- 81

648

m D

h

8 i 1 i élément cubique

(2.63)

Voyons quels éléments sont susceptibles de vérifier le principe d u m a x i m u m discret et sous quelles conditions. R e m a r q u o n s tout d'abord q u e les matrices de rigidité élémentaires vérifient toutes k.. = 0 ce qui est évident puisque

^

k

= zL n dN./dx dN Idx dx = D . dN./dx. d/dx ( Y N. ) dx or

C o m m e les matrices de m a s s e élémentaires sont à diagonale d o m i n a n t e stricte, la matrice globale M + a A t K sera à diagonale d o m i n a n t e stricte.

41

Enfin, puisque l'assemblage des matrices élémentaires n'entraînera que le recouvrement des termes a ^ et a^2 des matrices élémentaires, le principe du m a x i m u m discret sera vérifié si les termes non diagonaux de [k e

] + 9At [ m e

] sont < 0 .

(Nous ne nous intéressons ici qu'aux schémas implcites car M E F I S T O ne prend pas en compte à l'heure actuelle les avantages possibles d'une résolution explicite).

Il est évident que cette condition est toujours vérifiée pour l'élément linéaire avec la matrice [ m ^ ] (et c'est pourquoi l'exemple 2.3 ne présente pas d'oscillations) et qu'elle n'est jamais remplie avec les éléments d'ordre supérieur combinés avec la matrice masse diagonalisée.

Envisageons maintenant le cas où on utilise la matrice masse consistante.

La (les) condition(s) sont :

Elément linéaire

h

sou —

Elément quadratique

I _L_

D At

_ L _

4 10 9 /j2 10 9

< 0

19 39£>9Ai

h - < 0

1680 120 h

99 5 6 7 D e At

120 h

- 36 162J>6A¿

1680 ~ 120 h

En conclusion, le principe du m a x i m u m discret est vérifié, pour l'équation de la diffusion monodimensionnelle par les éléments suivants :

- élément linéaire et matrice masse diagonale V A í

42

- élément linéaire et matrice masse consistante si D A i / A " > 1 / 6 6

- élément quadratique et matrice masse consistante si

DM

1/40 6 < < 1/10

h

2

Remarque : le fait de ne pas vérifier le principe du m a x i m u m discret ne signifie pas que la solution oscillera systématiquement. Tout dépend en effet du second m e m b r e qui dépend lui-même des conditions initiales et aux limites.

Le traitement du problème de diffusion du § 2.3 avec des éléments quadratiques ou cubiques revient à imposer une solution initiale de la forme

T

II est donc raisonnable d'initialiser une solution correcte sur de tels

éléments si on veut limiter les risques d'oscillations. Pour l'élément quadratique, on peut par exemple démarrer avec la forme suivante :

1/4

Dans le cas bidimensionnel, on peut facilement montrer que seul le triangle linéaire vérifie toujours le principe du m a x i m u m discret à condition que tous les angles du maillage soient aigus. Il en résulte que ce principe sera toujours vérifié dans le cas non stationnaire à condition d'utiliser la matrice masse diagonale.

R e m a r q u e : le fait pour un élément de vérifier le principe du m a x i m u m discret ne saurait être considéré c o m m e une assurance sur la qualité de la solution, ainsi que le montre l'exemple 2.3 pour les temps courts et le schéma M D .

43

2.4.3 - CAS DE L'EQUATION DE DIFFUSION-CONFECTION

MONODIMENSIONNELLE

La matrice de convection élémentaire s'écrit (si U et D sont constants)

-

— 1

1

1

1 pour un élément linéaire

- 3

1

- 4

- 1

3

- 4

4 "

4

0

pour un élément quadratique

(on notera que la matrice de convection globale est antisymétrique).

O n montre alors facilement, à partir du théorème énonce en 2.4.2 que le principe du m a x i m u m discret est vérifié par l'élément linéaire si p - zJL <

2 e t n

u

'n

ne l'est pas par l'élément quadratique dans le cas stationnaire.

En suivant la m ê m e démarche, on trouvera des conditions reliant le nombre de Péclet numérique P n

et le nombre de courant U.At/h pour véfifier le principe du m a x i m u m discret dans le cas transitoire.

2.5 - C H O I X D U S C H E M A D E DISCRETISATION

2.S.1- LES CAUSES DES OSCILLATIONS NUMERIQUES

II existe trois causes principales des oscillations qu'on peut observer dans les solutions numériques de l'équation de diffusion-convection par la M E F G .

- U n e condition de Dirichlet sur un exhaure. L e remède consiste à changer cette condition en une condition de flux diffusif nul, ce qui est généralement plus réaliste. Dans le cas où une telle condition doit être imposée, l'utilisateur doit être conscient du fait qu'une couche limite est susceptible de se développer et prévoir en conséquence un maillage suffisamment fin pour "suivre" cette couche limite.

- U n nombre de Péclet numérique (P n

) trop élevé.

- D e s conditions initiales ou sur les limites amont très "raides" (par exemple du type échelon de température). La solution comporte alors trop de composantes spectrales dans le domaine des courtes longueur d'ondes qui ne peuvent être correctement suivies par un maillage trop grossier.

44

Dans les deux dernier cas, les oscillations sont produites par la dispersion numérique des schémas (i.e. les composantes spectrales de courte longueur d'onde ne sont pas propagées à la bonne vitesse). Le remède consiste alors à affiner la discrétisation spatiale. L'utilisation de schémas upwind supprimera ces oscillations mais au prix d'une diffusion numérique souvent très élevée et il n'y aura aucun avertissement (i.e. aucune oscillation) pour indiquer que la solution numérique est mauvaise.

Faut-il éliminer à tout prix les oscillations ? Il est certain qu'une solution numérique lisse est a priori plus satisfaisante mais elle est susceptible de coûter très cher dans le cas où les termes convectifs sont prépondérants si l'on veut éviter d'utiliser des schémas upwind. Si les oscillations apparaissent et persistent dans une zone considérée c o m m e importante du domaine d'étude, il est évident qu'il faut mailler plus finement cette zone. Si des oscillations dont l'amplitude reste faible n'apparaissent que dans des zones considérées c o m m e moins importantes, alors la solution qui n'oscille pas dans le reste du maillage présente probablement une bonne précision.

Enfin on notera que le schéma d'intégration temporelle est lui aussi susceptible d'introduire des oscillations (cf. 2.2 et 2.3).

2.5.2 - LE NOMBRE DE PECLETNUMERIQUE

Bien qu'il soit parfois vrai que le nombre de Péclet numérique ne doit pas dépasser 2 pour éviter les oscillations, c'est en général la combinaison d'un nombre de

Péclet élevé et de gradients importants dans la direction de l'écoulement qui produit les oscillations. Il est en fait admissible d'avoir des nombre de Péclet élevés dans certaines zones pourvu que les gradients y soient (ou y restent) faibles.

2.5.3 - LES SCHEMAS UPWIND

Ainsi que nous l'avons précisé à plusieurs reprises, les schémas up wind qui existent à l'heure actuelle doivent être maniés avec précaution. Les schémas up wind classiques consistent à décentrer progressivement la discrétisation spatiale des termes convectifs en fonction de l'importance du nombre de Péclet numérique, ce qui se traduit par l'introduction d'un terme de diffusion artificielle (numérique) permettant de limiter à 2 la valeur réelle de P n

. Ces schémas, qui ont le plus souvent

été développés pour traiter le cas stationnaire, s'avèrent donc trop diffusifs en régime transitoire où la limite (simpliste) de 2 pour la valeur de P n

est loin d'être imperative dans tous les cas.

45

La forme de l'équation de diffusion-convection utilisée dans la bibliothèque C V 2 D permet de générer un schéma upwind manuel en contrôlant la diffusion artificielle introduite. E n effet, le nombre de Péclet numérique dans la direction de l'écoulement est :

P =

d + r

IVI

II suffit donc d'augmenter la dispersivité longitudinale pour diminuer P n

et réduire les oscillations. Cette technique simple permet de ne pas introduire de diffusion transversale artificielle, ce qui est le cas de presque tous les schémas up wind qui existent à l'heure actuelle lorsque l'écoulement n'est pas parallèle aux frontières des éléments. Elle présente par ailleurs l'avantage de quantifier exactement la diffusion numérique qui a été ajoutée.

2.5.4 - DIAGONALISATION DE LA MATRICE MASSE

Nous avons montré précédemment que la réponse spectrale des schémas

M E F G utilisant une matrice de masse constante est bien meilleure que celle des m ê m e s schémas utilisant une matrice masse diagonalisée (ce qui explique que les résultats obtenus par la M E F G sont souvent bien meilleurs que ceux obtenus par les différences finies centrées pour un schéma d'intégration temporelle du type Exiler).

L'utilisation des schémas M D se justifie lorsque le schéma d'intégration temporelle est explicite. Dans la mesure où les blocs T R L C et T R L V du programme M E F I S T O ne sont pas conçus pour utiliser les avantages d'une résolution explicite, l'emploi de la matrice de masse consistante s'impose.

2.5.5 - CHOIX DES ELEMENTS

II est souhaitable, dans la mesure du possible, d'utiliser un maillage tel que les côtés des éléments sont parallèles ou perpendiculaires aux lignes de courant. C'est pourquoi il est préférable d'utiliser des quadrilatères.

La référence [5] analyse la réponse spectrale pour un problème de convection pure monodimensionnel des éléments linéaires et quadratiques et la compare aux schémas différences finies d'ordre 1 (upwind), 2, 4 et 6 (les différences finies centrées et upwind présentent la m ê m e réponse spectrale pour ce qui concerne la vitesse de phase mais le schéma upwind introduit un amortissement très important). Il apparaît que l'élément quadratique est le moins dispersif (il introduit cependant une onde supplémentaire qui se propage en sens inverse de l'écoulement

46

mais dont l'amplitude est généralement si faible qu'elle est difficile à détecter). Pour un problème bidimensionnel, l'élément quadratique Lagrange (type 5 de la bibliothèque C V 2 D ) est probablement le meilleur mais il est aussi plus sensible aux conditions initiales et aux limites.

2.6 - EXEMPLE

Nous nous proposons de visualiser sur un exemple simple les causes des oscillations. Pour ce faire, considérons le problème défini sur les figures 2.9 et 2.10.

Les lignes de courant ont été calculées à l'aide de la bibliothèque de diffusion D F 2 D

(fig. 2.13 et 2.14).

Nous utiliserons c o m m e maillage de base la grille représentée sur la figure 2.11. (On notera qua la vitesse n'est théoriquement pas définie sur les deux coins supérieurs de l'obstacle qui constituent des angles rentrant dans le domaine où prend place l'écoulement).

L'équation à résoudre est l'équation (1.1). La vitesse aux points d'intégration d'un élément est calculée à partir de la valeur de la fonction de courant aux noeuds. L e c h a m p de vitesse et le paramètre physique u) de l'équation (1.1) sont identiques dans tous les calculs qui suivent (w = 1).

La figure 2.15 montre les isothermes obtenues en régime stationnaire pour une condition aval de flux diffusif nul et pour les valeurs suivantes des paramètres physiques de l'équation 1.1

Avec ces valeurs, le nombre de Péclet numérique et suffisamment faible pour qu'il n'y ait aucune oscillation dans la solution calculée.

Diminuons maintenant, en gardant le m ê m e maillage, les dispersivités et le coefficient de diffusion :

Q

L

= 0 . 1 a r

= 0.02 d

0

= 0.01

Le nombre de Péclet est maintenant largement supérieur à 2 sur la majeure partie du domaine. O n observe alors des oscillations dans la zone située en amont de l'obstacle. Par contre, la zone située en aval de l'obstacle ne présente pas d'oscillations décelables (fig. 2.16).

47

Si nous modifions la condition sur la limite aval en imposant T = 0, des oscillations sévères apparaissent dans la zone située en aval de l'obstacle et la solution numérique est très mauvaise (fig. 2.17). En effet, le maillage est trop grossier pour suivre correctement la couche limite qui se développe au voisinage de l'exhaure.

Les figures 2.18 et 2.19 montrent les solutions obtenues dans les deux derniers cas avec un maillage beaucoup plus fin qui est représenté sur la figure 2.12.

Le nombre de Pédet numérique est partout inférieur à 2 et on observe aucune oscillation. En particulier, les couches limites sont parfaitement "suivies" par ce maillage.

La comparaison des solutions numériques obtenues avec le maillage grossier et le maillage fin dans le cas où la condition aval est 3 T / 9 n = 0 montre que la solution grossière est bonne dans la zone située en aval de l'obstacle. Dans la mesure où la solution en amont de l'obstacle nous intéresse peu dans le cas d'un problème de pollution, le maillage grossier s'avère suffisant pour traiter le problème

(fig. 2.20).

Les figures 2.21 a, b, c présentent la solution du régime transitoire avec les paramètres physiques utilisés dans le premier cas ((XL = 0.5, OLJ = 0.1, d o

= 0.1).

Bien que la solution en régime permanent n'oscille pas, la solution transitoire présente de légères oscillations pour les premiers pas de temps (non représentée sur les figures). Ces oscillations sont dues à l'échelon de température imposé sur l'obstacle qui fait intervenir trop de composantes spectrales dans le domaine des courtes longueurs d'onde. La vitesse de phase numérique de ces composantes est très différente de la vitesse réelle et la superposition des composantes de courtes longueurs d'onde produit les oscillations. Ces oscillations disparaissent cependant rapidement car le paramètre d'amortissement de ces composantes est élevé.

48

REFERENCES BIBLIOGRAPHIQUES

[1] BRGM - Programme MEFISTO

[Z] Donea (J.), Giuliani (S.), Laval (H.) 1979.- "Accurate explicit finite element schemes for convective-conductive heat transfer problem" in Finite

Element Methods for Convection dominated flows.- A S M E , A M D , vol. 34

[3] Gresho (P.), Lee (R.L.) 1979.- "Don't suppress the wiggles they're telling you something". In Finite Element Methods for Convection Dominated Flows.-

A S M E , A M D , vol. 34

[4] Varga (R.S.) 1962.- "Matrix Iterative Analysis" .- Prentice Hall

[5] Gresho (P.M.), Lee (R.L.), Sani (R.L.) 1978.- "Advection-Aduction dominated flows, with emphasis on the consequences of mass lumping".- In: Finite

Elements in Fluids, vol. 3, Wiley

6

è

Figure 2.1 - Solution de l'équation de diffusion-convection pour différentes valeurs de P E = U . L / D

o

L

\

O

«4-

Figure 2.2 - Facteur d'amplication (schéma d'Euler) o

51

Figure 2.3 - Rapport de la fréquence semi-discrète w à la frequence réelle w en fonction du nombre d'onde adimensionnel pour les éléments finis linéaires I D

Figure 2.4 - Rapport du paramètre d'amortissement semi-discrétisé au paramètre d'amortissement réel en fonction du nombre d'onde aHîmongirmn«»! pour les éléments linéaires I D

52 vu i.e 2.1 3.1

t.t

Figure 2.5 - Réponse fréquentielle des schémas d'Exiler implicites convection pure)

H 1 1 h

1 1 1 1

1 1 1 i.e 2.1

1 1

l.t S.t

Figure 2.6 - Amortissement introduit par les schémas d'Exiler implicites

(convection pure)

53 t- 0.4 r

1O

U

Figure 2.7 - Evolution temporelle des températures nodales pour l'exemple 2.3

D.1

i

0.8

0.7

0.9 m

-i

1

I

1

1

II«

M

1

I'

\

\\

l\

>(

Vk¿««» n o _

0.6

0 . 5 r-

0.4

0.3

0.2

0.1

I

•1

i

1

\

1 \

\

0.004

Í

-

—

\ y

i X\ i i !

i \

A

V

*\

x\

^ t = 0.04

OS ^ ^ t L »

-

-

—

Figure 2.8 - Profils de température pour l'exemple 2.3

1

0.2

1 1 , 1 ,

0.4

0.6

0.B

1.0

Figure 2.9 - Conditions aux limites pour le problème hydraulique

T Ï O

•fcT/fcnïO

T ï . 4

Tt -t Til

Figure 2.10 - Conditions aux limites pour le problème thermique eu.

T : O

2 . 0 -

- 3 . 0 - 2 . 0 - 1 . 0 0 . 0 1.0 2.0 3.0 4 . 0 5.B

- 3 . 0 - 2 . 0 - 1 . 0 0 . 0

Figure 2.11 - Maillage grossier

55

2 . 0 -

- 3 . 0 - 2 . 0 - 1 . 0 B . 0 lTe 2.0 3.B ¿.0 5.0

1-

-

-

—

1

_

-

7.0

1.0l.B

- 3 . 0 - 2 . 0 - 1 . 0 0 . 0 1.0 2.0 3.0 4.0 5.0

-0.0

Figure 2.12 - Maillage fin

Figure 2.13 - Lignes de courant calculées avec le maillage fin

(de 0.1 à 0.9 par pas de 0.1)

Figure 2.14 - Lignes de courant calculées avec le maillage grossier de 0.1 à 0.9 par pas de 0.1)

56

Figure 2.15 - Solution stationnaire - maillage grossier d o

= 0.1,

O L

= 0.5, ctj = 0.1

isothermes de 0.1 à 0.9 par pas de 0.1

i :

^ ^ —

Figure 2.16 - Solution stationnaire -maillage grossier do = 0.01, CXL = 0.1,

O T

= 0.02

isothermes de 0.1 à 0.9 par pas de 0.1

Figure 2.17 - Solution stationnaire - maillage grossier do = 0.01, a

L

= 0.1, ct

T

= 0.02

isothermes de 0.1 à 0.9 par pas de 0.1

TrO

57

Figure 2.18 - Solution stationna ire - maillage fin d o

= 0.01,

O L

= 0.1, a-p = 0.02

isothermes de 0.1 á 0.9 par pas de 0.1

£3

.«O

DP

Figure 2.19 - Solution stationnaire -maillage fin do = 0.01,

OIL

= 0.1,

<*T

= 0.02

isothermes de 0.1 à 0.9 par pas de 0.1

Figure 2.20 - Comparaison des solutions stationnaires obtenues avec le maillage grossier ( ) et le maillage fin ( )

58

o.-i

Figure 2.2 l.a - Solution transitoire - maillage grossier do = 0.1,

(XL

= 0.5, aj = 0.1 ; t = 0.05

isothermes de 0.1 à 0.9 par pas de 0.1

Figure 2.21.b - t = 0.55

Figure 2.21.C - t = 3.05

59

là

•

Figure 2.2l.d - t = 8.05

Figure 2.21.e - t = 18.05

réalisation service reprographie du B R G M

86 SGN 194 EA