application de la méthode des éléments finis à la - Infoterre
Contents
1. eel AB Figure 2 4 5 scheme MC pecheme MD 4 by Riz Rapport de la fr quence semi discr te w la fr quence r elle w en fonction du nombre d onde adimensionnel pour les l ments finis lin aires 1D th mes nc 1 2 bcheme MD 8 19 0 2 6 3 8 40 0 58 6 68 8 78 D 68 8 58 1 8 RER MERS EEE Rh ak Pj Rih wR Rapport du param tre d amortissement semi discr tis au param tre d amortissement r el en fonction du nombre d onde adimensionnel pour les l ments lin aires 1D eller El st 32 t 2 8 8 z e 0 S 0 08 6 04 4 0 r e4 a w At 1 8 2 8 3 2 4 8 5 8 8 0 Figure 2 5 R ponse fr quentielle des sch mas d Euler implicites convection pure 8 28 wo DE p 1 8 2 8 3 8 4 8 5 8 B B Figure 2 6 Amortissement introduit par les sch mas d Euler implicites convection pure 53 1 0 gt autren arnete we recreme TIC a2 re 1074 1073 1072 107 10 t Figure 2 7 Evolution temporelle des temp ratures nodales pour l exemple 2 3 me m oRaction ar nacte acre me WC 0 7 ne ES TS mno Figure 2 8 Profils de temp rature pour l exemple 2 3 54 Figure 2 9 Conditions aux limites pour le probl me hydraulique D Tl n 0O Figure 2 10 Conditions aux limites pour le probl me thermique 3 0 2 1 8 8 2 1 2 2 8 3 2 4 0 S B 2 0 2 8 ue Re D SOUT Rat S E V EE E MU OE SAE NORTE 1 2 ere ee E E al Gaa2. 2 OF AI a 2 SJ taz ay la 2 on on on dy dy 1 32 ot J est la matrice jacobienne de la transformation g om trique Les termes de J sont obtenus en d rivant par rapport et n la relation 1 31 On obtient par exemple ox oN AR en CEA dE l Ox aN I ar an i 1 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 une int gration plus simple sur l l ment de r f rence VT f x dx dy f x E n det J dE dn 1 33 det J tant le d terminant de la matrice jacobienne J 1 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 ls f E d Edn Z ner o amp 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 On notera que les tranformations effectu es et la technique d int gration utilis e n cessite le calcul des expressions du type Nj E ANi E 8E Ces expressions ne d pendent pas de la g om trie de l l ment r el et peuvent donc tre
3. o P U L D est le nombre de P clet qui mesure l importance du transport convectif par rapport au transport diffusif Tant que Pe est faible la solution varie doucement sur le domaine mais lorsque P croit la solution devient du type couche limite c est dire T x To sauf au voisinage de l extr mit x L o elle s adapte brutalement la condition d exhaure T L Ty fig 2 1 L paisseur de cette couche limite peut tre estim e en galant les termes convectifs et diffusifs de l quation sur une faible distance 6 soit U AT DAT 5 En fin de texte 24 ce qui conduit 6 L 1 P qui indique que l paisseur de la couche limite est approximativement inversement proportionnelle au nombre de P clet La difficult de ce probleme r side dans le fait que toute la physique de la solution est contenue a l int rieur de cette mince couche limite lorsque P gt gt 1 La discr tisation de l quation 2 1 par la MDFC ou la MEFG lin aire i e avec des l ments lin aires sur une grille uniforme de pas h donne au noeud i es a Ls i 1 T _ 2T T D cae EAA soit 2h h Pi eg e Ae E 2 5 n h P rires 2 6 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 Z5 Lar 2 Pn in es FES K 2 7 i ee La forme de ce
4. pour r soudre l quation de diffusion convection La 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 DF2D cf 1 le champ de vitesse calcul est continu sur l int rieur des l ments et discontinu la travers e des fronti res On notera cependant que le calcul de la matrice de rigidit est r alis en 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 a 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 sym trique s effectue de la mani re suivante calcul pr alable des det J lt 8N 3x gt lt 8N
5. quations d coupl es AE e ar 2 18 On peut alors examiner l influence de en tudiant l quation simple 29 VAV S0 2 19 ce qui permet de comprendre le comportement des diff rents modes qui forment la solution du probl me initial La solution analytique de 2 19 est i 1 _ l Vit exp A t V erp AAt V 2 20 L quation 2 19 discr tis e par le sch ma d Euler donne 1 68AA V l 1 1 0 AA1 V soit vittca y 2 21 ou ee Gs ee me er 2 22 p 1 OA est appel facteur d amplification Dans le cas o le sch ma d int gration temporelle fait intervenir la valeur de la variable plus de deux instants A est une matrice matrice d amplification L tude de la stabilit de la dispersion du sch ma est alors bas e sur les valeurs propres de la matrice d amplification La solution exacte 2 20 indique que VE perv sf SO vitt vi gi A 0 2 23 compte tenu des caract ristiques du syst me 2 13 les modes sont r els et 2 0 Le sch ma sera donc stable si A lt 1 la deuxi me condition est satisfaite si 0 soit 1 1 AAe Se ET ET i 2 24 L in galit de droite est toujours satisfaite pour toutes les valeurs possibles des param tres At gt 0 et 0s 6 lt 1 et l in galit de gauche est v rifi e si 9 2 0 5 Si 6 lt 0 5 l in galit de gauche impose 30 AAt lt 2 1 286 8 lt 1 2 2 25 ce qui pour donn fixe une val
6. 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 CV2D Soit e un l ment de la triangulation poss dant n noeuds La solution approch e sur cet l ment s crit M a c a un 1 28 l ka li Les coefficients Cj 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 variationnelle 1 24 s crit sur cet l ment VaN dv WaC 80 N dV N hdS L y I L S q VN av I V V Compte tenu de la d finition des flux 1 7 et 1 8 et de 1 28 l quation pr c dente donne UNN av C VN D YN ava N vUN pav c QN dV e io J J e 1 j e i j J ye V V V 1 29 h N dS S Cette quation peut s crire sous la forme matricielle suivante m iC PH TRI 2 iC path b th 1 30 avec Cj valeurs nodales de la concentration G diC dt m ye wN N av _ matrice de masse k le VN DYF N j dV matrice de diffusion dispersion z N Y UN dV matrice de convection fy QN dV vecteur des sollicitations r parties h N dV ive eo f T vecteur des flux diffusifs impos s sur q i les fron
7. celle de l inconnue variationnelle c est dire la temp rature ou la concentration La vitesse en un point de l l ment est alors calcul e partir de ses valeurs aux noeuds de l l ment par U U N x XJ J U ti y WY J Le champ de vitesse ainsi d fini est continu sur le domaine d tude V I n est cependant qu une approximation du champ 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 Ce champ 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 polynomes 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
8. 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 L F gt Osi F gt 0 40 Les l ments d ordre sup rieur ne v rifient pas le principe du maximum qe j 6 eee T amp 28 en Se Le ag pour pean Anagr vane 8 8 16 me a a a 2 59 h pour l l ment cubique et L i _ D m _ 567 162 2 60 120 h 444 162 367 S 1296 891 AEE EA nas ym 1296 La pr sence de termes non diagonaux gt 0 dans ces matrices l mentaires va entra ner la pr sence de termes non diagonaux gt 0 dans la matrice globale Dans le cas transitoire le principe du maximum discret sera v rifi si la matrice M aAtK v rifie le th or me pr c demment nonc Les matrices de masse l mentaires consistantes et diagonalis es des l ments lin aires quadratiques et cubiques sont hIL 1 h 1 0 m m l ment lin aire 2 61 CT 611 2 D 240 A eee 4 1 2 e _ h 1 eC 30 sym 2 D Te ei l ment quadratique 2 62 16 4 h 128 19 99 36 e _ h 1 mo 1680 128 36 99 MDI g 1 l ment cubique 2 63 sym 648 8l 3 648 3 Voyons quels l ments sont susceptibles de v rifier le principe du maximum discret et sous quelles conditions Remarquons tout d abord que les matrices de rigidit l mentaires v rifient toutes gt k 0 ce qui
9. de l obstacle et la solution num rique est tres 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 clet 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 9T dn 0 montre que la solution grossi re est bonne dans la zone situ e en aval de l obstacle Dans la mesure ou 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 a 0 5 ay 0 1 d 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
10. el fait lui m me intervenir le calcul du jacobien au point d int gration L enchainement des op rations est le suivant Pour chaque point d int gration calcul de la matrice jacobienne 2 x 2 partir du tableau VNI et du tableau VCORE 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 JACOB calcul si n cessaire des quantit s A8N 9x et 9Nj dy au point d int gration On d duit de la relation 1 32 oN aN Or aE a sii La 1 34 y on Cette op ration est r alis e par le sous programme DNIDX partir du tableau VNI et de l inverse de la matrice jacobienne qui vient d tre calcul e Le r sultat est stock dans le tableau VNIX qui contient s quentiellement lt aN ax gt lt aN ay gt pour E E lt aN dx gt lt aN dy gt pour f 6 1 34 soit n x 2 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 VNI et ou VNIX ainsi que les tableaux contenant les constantes de milieu et de force de l l ment Cette organisation a t retenue car les tableaux a stocker du type VNIX 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
11. est vident Lj puisque J gt k FP ON ax aN ax dx D ON 0x a ax oy N dx J J J or gt N J Comme les matrices de masse l mentaires sont diagonale dominante stricte la matrice globale M aAtK sera a diagonale dominante stricte 4 Enfin puisque l assemblage des matrices l mentaires n entra nera que le recouvrement des termes aj et a22 des matrices l mentaires le principe du maximum discret sera v rifi si les termes non diagonaux de k 6At m sont lt 0 Nous ne nous int ressons ici qu aux sch mas implcites car MEFISTO ne prend pas en compte l heure actuelle les avantages possibles d une r solution explicite I est vident que cette condition est toujours v rifi e pour l l ment lin aire avec la matrice mP 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 D At 1 k Pen lt 0 soit 6 El ment quadratique h D At Y j D At 1 _ lt 0 sit lt lt 30 3h 4 100 h2 108 2h B8DeAt a lt 30 3h El ment cubique 19 39D0At 1680 120k 99 567 DO At _ o ee EE m 1 a pas olution 1680 120 RS er 36 162D04A1 lt 1680 120h En conclusio
12. 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 g U C U C n u ay ce qui tend imposer C 0 Ce 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 2 1 Par contre une condition en flux diffusif nul l exhaure implique si D est isotrope dC dn 03 ce qui appara t coh rent avec l hypoth se que l exhaure se trouve dans une r gion ou les gradients sont tr s faibles 1 4 DISCRETISATION SPATIALE ET FORMULATION MATRICIELLE Nous utiliserons la formulation en flux diffusif associ e a 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 somme 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
13. 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 d e I appara t finalement que le choix d une formulation variationnelle 1 20 0 sur S 1 27 int ou 1 24 conduit dans le cadre des fonctions de classe C9 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 a travers les fronti res des l ments Dans la formulation en flux diffusif les conditions sont quivalentes mais s appliquent au flux diffusif 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 EAU 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
14. la solution ainsi que le montre l exemple 2 3 pour les temps courts et le sch ma MD 43 2 4 3 CAS DE L EQUATION DE DIFFUSION CONVECTION MONODIMENSIONNELLE La matrice de convection l mentaire s crit si U et D sont constants pour un l ment lin aire 4 er pour un l ment quadratique on notera que la matrice de convection globale est antisym trique On montre alors facilement partir du th or me nonce en 2 4 2 que le principe du maximum discret est v rifi par l l ment lin aire si p Ur lt 2 et qu il En D ne l est pas par l l ment quadratique dans le cas stationnaire En suivant la m me d marche on trouvera des conditions reliant le nombre de P clet num rique P et le nombre de courant U At h pour v fifier le principe du maximum discret dans le cas transitoire 2 5 CHOIX DU SCHEMA DE DISCRETISATION 2 5 1 LES CAUSES DES OSCILLATIONS NUMERIQUES Il existe trois causes principales des oscillations qu on peut observer dans les solutions num riques de l quation de diffusion convection par la MEFG Une condition de Dirichlet sur un exhaure Le 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
15. m me boucle sur les l ments la matrice de rigidit puis la matrice de masse C est pourquoi les 10 premi res variables du common LMTRAV sont des indicateurs qui sont mis z ro par le programme MEFISTO 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 1 5 5 CALCUL DE LA MATRICE DE RIGIDITE L quation effectivement trait e dans la biblioth que CV2D suppose que la divergence de la vitesse est nulle sur le domaine soit VU 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 do et les dispersivit s ar et ar sont constantes sur l l ment mais que la vitesse est variable Suivant les conventions adopt es dans MEFISTO et d taill es dans 1 do ap et ar sont des constantes de milieu la vitesse est une propri t nodale dont il faut sp cifier l interpolation sur l l ment Un choix naturel consiste adopter une interpolation de la vitesse identique
16. 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 De 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 Il 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 i
17. 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 MEFISTO 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 mode 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 EAU Les principales tudes r alis es l aide de MEFISTO 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 tro combustion les transferts de substance dissoute en aval d un stockage souterrain de gaz naturel pompage par un puits incomplet SOMMAIRE INTRODUCTION 1 LA BIBLIOTHEQUE D ELEMENTS CV2D 1 1 L quation de diffusion convection 1 2 Notations et conditions aux limites 1 3 Formulation variationnelle 1 3 1 Formulation en flux total 1 3 2 Formulation en flux diffusif 1 4 Discr tisation spatiale et formulation matricielle 1 5 Calcul effectif des matrices l mentaires 1 5 1 Transformation des d riv es premi res 1 5 2 Transformation des int grales 1 5 3 Int gration num rique 1 5 4 Organisation du calcul 1 5 5 Caleul de la matrice de rigidit 1 5 6 Calcul de la matrice de m
18. t s positifs et n gatifs Le saut de qp au point x est d fini par l l I Il est clair que le saut qp 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 i e ind finiment continues et d rivables En appliquant le th or me d Ostrogradski l quation 1 5 on obtient q W av pwav q n WdsS 1 19 S V U D composons l int grale sur la fronti re du domaine qn was q was q WdS S l S n S n c q Puisque W est nulle sur So et compte tenu de 1 13 ce terme s crit finalement q n WdS W l dS S S q et l quation variationnelle quivalente 1 15 est P a W av pwav W IdS 1 20 Ns s q Dans la m thode des l ments finis les int grales portant sur le domaine sont remplac es par la somme des int grales sur chaque l ment de la triangulation Si les fonctions W et C sont de classe C l quation variationnelle 1 20 est quivalente W a Nav Wq l dS 0 s 1 21 Comme 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 CO et plus exactement telles que leurs d riv es soient des continues sur les fronti res internes de l
19. tude et qu il n y ait pas de source ou de puits de cette substance Le principe du maximum nonce alors que la concentration maximum 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 maximum 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 aT ae ie a sur V aT an 0 sur Se T g sur S 2 56 T x 0 E sur V v rifient aussi le principe du maximum qui s crit max max T max g t max 0 max Q Mie i 7 2 57 t min 0 min_ Q lt T x 0 lt V min mina min g avec E Sq x 0 T i 2 4 2 CAS DE L EQUATION DE LA DIFFUSION On peut montrer facilement que l l ment lin aire monodimensionnel v rifie le principe du maximum 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 ajj est diagonale dominante avec ajj lt 0siifjet ajj gt Osiie 1 n alors A7 gt 0 une matrice n x nest dite diagonale dominante si X a 20 Vi J La matrice de rigidit l mentaire de l l ment lin aire monodimensionnel pour l quation de la diffusion s crit l 1 2 58 o h est la longueur de l l ment et D le coefficient de diffusion Il est
20. 9y gt aux points de Gauss si n cessaire cf 1 5 4 r 0 Boucle sur les points d int gration Calcul de la vitesse Ux Uy au point d int gration de coordonn es amp n 3options sont possibles a on conna t la vitesse aux noeuds de l l ment alors U bs Ti N amp n J 1 34 U gt T N 8 n J b on conna t le potentiel hydraulique aux noeuds et le tenseur de perm abilit constant sur l l ment On suppose que l interpolation du potentiel h sur l l ment est identique celle de l inconnue variationnelle soit h gt h N ce qui est le cas si le calcul de h a t effectu avec la biblioth que DF2D On a alors d apr s la loi de Darcy U H lt Kh aN Jaxt K x Xx p ON fay U K y xy N A J h aN ox K gt h N y J J yy J J o Kyy Kxy et Kyy sont les composantes distinctes du tenseur de perm abilit sur l l ment c On conna t la fonction de courant On suppose alors que y D y N J sur l l ment et que la vitesse est U apla gt y aN Jay 5 U aw ax Ss w aN Jax j en plan U 2 aw ay en radial o r est le rayon r x ou y selon la s valeur de NAXIS dans le common PROB U aw ax cf 1 Calcul du tenseur D au point d int gration l aide de 1 2 Modification de r pour i l netj 1 n r r aN x D 8N x D N y f U N LJ L J XX l XY l x i aN 3y
21. D 3N 3x D N l f U_ N J xy t YY E Y J L det J E n Pk o Pk est le poids du point d int gration 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 MILI 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 Ce 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 D 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 m F wN N aV w 2 N Eon N Eng de J en P 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 NMDIAG du common PROB en appelant le bloc PROB du programme MEFISTO La matrice de ma
22. X IU O ao aU Ex DR d vy IUI o 1 2 a a U_ U xy ye UI UI UT UT 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 Ya T N UTE Y dt lt A Ye avec T temp rature Qr terme source U vitesse de Darcy de l coulement Ya capacit calorifique du milieu tenant compte de la pr sence de la phase fluide f _ capacit calorifique de la phase fluide 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 a 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 sont Ap A B yp 1 4 oe A A Br Vp o Br et Br sont respectivement les dispersivites thermiques longitudinales et transversales L quation 1 3 est formellement identique l quation 1 1 condition d identifier w et Ya Yf do et o Yf GL et aT BL et BT Q et Qr Yys 1 2 NOTATIONS ET CONDITIONS AUX LIMITES La biblioth que d l ments CV2D 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 nomb
23. a triangulation i e sur Sint l quation variationnelle 1 20 est quivalente Wia f dv W Das W dS 0 gt As f a S la d 1 22 e h int L quation 1 22 montre que dans la classe des fonctions CO 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 Iq 0 sur os 1 23 1 3 2 FORMULATION EN FLUX DIFFUSIF Dans cette formulation le th or me d Ostrogradski n est appliqu qu au terme de flux diffusif En r crivant 1 15 sous la forme a q f WdV 0 U et en tenant compte de gf WdV gt W dV gi Wn dS Pa L L L S L l L U on obtient awidv af wav WfdV gin WdS vo as v StS an soit en tenant compte de 1 14 et puisque W est nulle sur S qo W dV qi W dV WidV Wh dS 1 24 U i S U V q Si les fonctions W et C sont de classe C l quation variationnelle 1 24 est quivalente Wig PdV W g h dS 0 1 25 U J S h et les quations d Euler Lagrange sont 1 10 et 1 14 Si W et C ne sont que de classe CO l quation 1 24 conduit pa e h int Les quations d Euler Lagrange sont alors wa nav Wal has Wig dS 0 1 26 i S la restriction de l quation 1 10 l
24. aires 4 noeuds et biquadratiques incomplets 8 noeuds cf 2 2 3 OSCILLATIONS AVEC L EQUATION DE 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 interpret s comme un avertissement 36 Consid rons le domaine 0 lt x lt 1 discr tis sur une grille 4 pas constant h form e d l ment lin aires et le probl me aT _a T 2 46 ot 7 se T 0 t 0 T l x 0 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 07 Pour l quation de la diffusion et un l ment gt wr h 2 1 m a 1 9 dans le cas MC matrice masse consistante 2 47 h 1 0 m aloa dans le cas MD matrice masse diagonale 2 48 ne D 1 1 o D est le coefficient de diffusion 2 49 oh 1 1 D 1 dans le probl me trait 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 soit 4 1 T 22 T l 1 T 1 T 0 2 ped A Te 2 50 Ty Ty dans le cas d un sch ma MC o T est
25. ant 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 CV2D7Q de la biblioth que CV2D dans lequel l int gration est r alis e avec 3 points Dans le cas axisym trique on a dS 2fr dr dz et l int grale valuer est 2n RON E r CE Je CE d 1 Nous donnons ci dessous l expression du vecteur f aux noeuds d une ar te droite Ar te lin aire le flux varie lin airement f L 2 1 hy se liz la empan E L 3 ri To rot To h Ten 12 h cas axisym trique o rj est la valeur du rayon aux noeuds Ar te quadratique droite h3 2 attention h2 est la valeur du flux au noeud 3 de l ar te 3 h2 Le RUN le flux varie de fa on quadratique hase kook o les N k sont les fonctions h j oe d interpolation de l l ment mono dimensionnei quadratique f Z L 4 1 2 hy s 30 1 4 2 h 3 2 2 16 h en plan 2 pac Tritro ir tr 2r h s 30 r tra r 7r 27 h _ D i 3 f 3 en g om trie axisym trique 27 279 16 r ro he 21 1 6 UTILISATION DE LA BIBLIOTHEQUE CV2D La structure de la biblioth sque CV2D 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 ELEMLB appel des di
26. application de la m thode des l ments finis la simulation des transferts dans les eaux souterraines BRGM notice et mode d emploi du programme MEFISTO STEF ml 3 ll ro application de la m thode des l ments finis la simulation des transferts dans les eaux souterraines notice et mode d empioi du programme MEFISTO STEF M Recan octobre 1986 86 SGN 194 EAU BUREAU DE RECHERCHES GEOLOGIQUES ET MINIERES SERVICE GEOLOGIQUE NATIONAL D partement Eau B P 6009 45060 ORLEANS CEDEX 2 T l 38 64 34 34 RESUME Le pr sent rapport comprend la description le mode d emploi et un jeu d essai du programme MEFISTO STEF 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 MEFISTO et MEFISTO STEF 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
27. asse 1 5 7 Calcul des sollicitations r parties 1 5 8 Calcul des sollicitations sur les fronti res 1 6 Utilisation de la biblioth que CV2D NOTIONS SUR LES SCHEMAS UPWIND LA STABILITE ET LA DISPERSION DES SCHEMAS NUMERIQUES 2 1 Les sch mas upwind 2 1 1 Cas monodimensionnel 2 1 2 Extension 2 ou 3 dimensions Page 10 12 12 13 13 13 15 18 19 19 21 23 23 23 26 Page 2 2 Notions de stabilit dispersion et dissipation des sch mas num riques 27 2 2 1 Stabilit quation de la diffusion 21 2 2 2 Dissipation dispersion 31 2 3 Oscillations avec l quation de la diffusion 35 2 4 Le principe du maximum 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 Le nombre de P clet 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 De nombreux probl mes physiques sont mod lis s par l quation de diffusion convection Mis a 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 d
28. calcul es une fois pour toute pour chaque type d l ment 1 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 GAUSS stockage dans les tableaux VCPG et VKPG des poids et des coordonn es des points d int gration sur l l ment de r f rence carr GAUSST idem pour l l ment de r f rence triangulaire FONC NI 2D stockage dans le tableau VNI des valeurs Ny Ep SNi Er 9E INiEr n Le tableau VNI contient s quentiellement lt N gt lt 0N 66 gt lt 0N an gt pour lt N gt lt dN 06 gt lt N n gt pour b 1 34 lt N gt lt ON 06 gt lt 4N dn gt pour avec lt N gt 8 lt N 8 N 8 N 8 gt lt dN 1df gt 8 lt 9N 8 108 0N 6 98 gt 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 CV2D les quantit s calculer sont det J E qui appara t lors de la transformation des int grales et JN E 9 Le Le calcul des d riv es des fonctions d interpolation sur l l ment r
29. e 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 Ce rapport pr sente la biblioth que d l ments CV2D bas e sur la formulation de Galerkine Coupl e au programme MEFISTO 1 cette biblioth que permet de traiter un probl me de diffusion convection La formulation utilis e et le mode d emploi de la biblioth que CV2D 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 apparaitre des probl mes similaires a 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 met en vidence le
30. e la solution pour n importe quelle valeur du nombre de P clet Le rem de aux oscillations consiste suivre correctement la couche limite lorsque P 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 De 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 I 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 Il 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 clet lev appara t aussi dans les cas 2D ou 3D 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 La v r
31. e masse consistante s impose 2 5 5 CHOIX DES ELEMENTS Il 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 me r ponse spectrale pour ce qui concerne la vitesse de phase mais le sch ma upwind introduit un amortissement tr s important I 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 CV2D 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 DF2D fig 2 13 et 2 14 Nous utilise
32. e pour des valeurs sup rieures le signe de Al change a chaque pas Pour les algorithmes inconditionnellement stables 2 1 2 la valeur asymptotique du facteur d amplification A v rifie A o lt 1 Si 6 gt 1 2 alors quel que soit AAt Al lt 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 8 1 2 ou est tr s pr s de 1 2 et AAt gt gt 1 alors A 1 et les composantes modales lev es se comportent comme 1 i Ces oscillations dans le temps se manifestent souvent dans les calculs On peut dans ce cas filtrer ces composantes lev es en utilisant la moyenne de deux instants Vi 1 Vi 2 puisque 1 1 0 Remarque Si a 1 et At gt 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 tl par leur valeur exacte dans l quation discr tis e v rifie e t lt C Atk l pour tout t C est une constante ind pendante de At et k gt 0 est l ordre de pr cision ou taux de convergence On montre par ailleurs facilement que le sch ma d Euler v rifie k 1 pour tout O 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
33. e w est w sin kh kh w 1 2r coskh 1 2 38 et que le rapport du param tre d amortissement semi discr tis t au param tre d amortissement exact amp est Z 2 1 cos kh kh C 2r coskh 1 2 39 Le rapport des fr quences 2 38 et le rapport d amortissement 2 39 sont fonctions du nombre d onde adimensionnel kh Le r sultat important qui appara t sur la figure 2 3 est que les deux discr tisations MD et MC r duisent les fr quences Le sch ma MD exhibe une r ponse en f r 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 MC l augmente alors que le sch ma MD 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 sch ma d Euler l guation 2 33 o B i w Z ce qui donne o t dt 1 1 8 At iw amp SS SR 2 40 p t 1 6At iw amp Par analogie avec la solution exacte de 2 33 qui est donn e par 2 37 nous re crivons o t dt en iw T Solt ce qui donne compte tenu de 2 40 EE A E 1 a At iw 1 6At iw OQ dont on d duit que ig D At b a 2 41 35 et e l At Va b c 2 42 avec 1 1 AI 6 At 6 1 8 w At 2 43 w AT 2 44 2 45 c 14 6Z At 07 w At Le rapport w w dans le cas de la convectio
34. es 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 52 totale qui doit tre consid r e Pour fixer les id es nous consid rons l quation de diffusion convection monodimensionnelie homog ne sans terme source 2 T T ne d ox ay ou U et D sont constants Ce syst me est non dispersif Pour une solution initiale de la forme T 0 LT LRx x 0 e 2 29 la fr quence associ e est w kU et le param tre d amortissement est Dk 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 l U D a paces Se eme es 0 Bag OE PE a oy ea Psy cee 2T T 4 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 somme des N m E N N av N EN av N av l termes de son rang 5 m l quation pr c dente s 1J l quation pr c dente s crit r T D r 2T T J 2h ae ae Red 2 Pad a J 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 U D CR aN ere a ear dl 2 30 ou L est l op rateur d fini par LT Tj 1 Tj Tj 1s r 0 pour la cas MD matrice ma
35. eur maximale du pas de temps Le sch ma d Euler est dit conditionnellement stable pour 8 lt 1 2 et inconditionnellement stable pour 9 gt 1 2 Remarque Il est imp ratif de v rifier la condition A lt 1 En effet la solution de 2 21 peut s crire V Al VO et si A est gt 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 lt 2 1 28 doit tre v rifi e pour tous les modes i e tous les d k 1 n du syst me 2 13 Le plus grand mode max impose la restriction la plus s v re sur le pas de temps maximal En fait dans le cas de l quation de la diffusion on peut montrer que max O h7 ie est de l ordre de h o h mesure le pas de la discr tisation spatiale Le pas de temps maximum doit donc v rifier At lt constante he 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 Remarque Si le syst me 2 13 est issu d une discr tisation par l ment finis alors Amax est inf rieur la valeur propre maximale des l ments pris individuellement Le comportement de A en fonction de AAt est repr sent sur la figure 2 2 pour diff rentes valeurs de 6 La valeur de AAt pour laquelle A 0 est appel e la limite d oscillation parce qu
36. ff rents types d l ment ELEMOI CV2D ex cution des fonctions l mentaires pour l l ment de type 1 triangle lin aire 3 noeuds ELEMO2 CVZD ex cution des fonctions l mentaires pour l l ment de type 2 quadrangle lin aire 4 noeuds ELEMO03 CV2D ex cution des fonctions l mentaires pour l l ment de type 3 triangle quadratique 6 noeuds ELEM04 CV2D ex cution des fonctions l mentaires pour l l ment de type 4 quadrangle quadratique incomplet 8 noeuds ELEMOS5 CV2D ex cution des fonctions l mentaires pour l l ment de type 5 quadrangle quadratique Lagrange 9 noeuds CV2D3 calcul de la matrice de rigidit pour un des l ments de la biblioth que CV2D5 calcul de la matrice de masse CV2D6 calcul du r sidu r C CV2D7 calcul du second membre Q N dV CV2D7L calcul de h N d Sur ite arate lin aire CV2D7Q calcul de gt sur une ar te quadratique CV2D8 calcul et impression les gradients aux points d int gration CV2D9 calcule les vitesses aux points de Gauss et les stocks dans le tableau des propri t s l mentaires pour utilisation future par CV2D3 et CV2D6 CV2DNX 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 LMCV2D OLB Les sous programmes g n raux de MEFISTO se trouvent dans la biblioth que MEF OLB Le programme principal qui r a
37. i 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 a Wav q W av q Wn dS 1 17 i S U 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 diffusif 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 Avant d expliciter les deux formulations nous introduisons la notation suivante Soit x un point de Sint C est dire situ sur une fronti re interne de la triangulation On appelle de fa on arbitraire un c t de Sint le c t positif et l autre le c t n gatif Soient n et n les normales unitaires S au point x qui sont dirig es vers les directions positives et n gatives Il est vident que nt n Soient q et q les valeurs de qj obtenues en approchant du point x par les c
38. itable 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 clet 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 On 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 tres efficaces dans certaines situations 2 2 NOTIONS DE STABILITE DISPERSION ET DISSIPATION DES SCHEMAS NUMERIQUES Les sch mas upwind bri vement pr sent s dans le
39. 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 TT coth Pn 2 RP dk con hm PAC ET 2 11 ou Pj U h D est le nombre de P clet de maille sur l l ment i L examen de la relation 2 10 montre que Pp Pp lorsque Ph lt lt l et Ph 2 lorsque P gt gt 1 Par ailleurs la diffusivit originelle D de l quation 2 1 est remplac e par la diffusivit effective _ Pn Pn D D coth 9 2 2 12 qui donne D gt D lorsque P lt lt 1 et D gt 1 2 U h lorsque Pn gt gt 1 On voit que lorsque P est grand P gt 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 a 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 gt 2 peuvent tre interpr t es comme 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 amen croire la validit d
40. la temp rature du ler d l libre i e x 1 h Dans le cas du sch ma MD on a M hl ou I est la matrice identit la discr tisation MD est dans ce cas strictmeent identique au sch ma diff rences finies pusique les seules conditions aux limites sont du type Dirichlet On peut montrer cf 3 pour les r f rences que l inverse de la matrice masse consistante poss de les propri t s suivantes 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 La valeur de T t 0 est donn e par M lf Comme seul le premier l ment de f est non nul T t 0 est gal 1 h fois la premi re colonne de M71 soit T t 0 Mij h Compte tenu des propri t s de M7 rappel s ci dessus T 0 est gt 0 si j est impair et lt 0 si j est pair l amplitude de T 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 Dans le cas du sch ma MD matrice masse diagonalis e on a T 0 1 h2 et T dans la mauvaise direction Il semble donc premi re vue que la m thode de 0 0 j 2 N c est dire qu aucune des temp ratures nodales ne d marre Galerkine consistante pr te caution Comparons la s
41. leur moyenne de pet pnest une bien meilleure approximation de n f 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 On observe ce qui a t pr vu pr c demment savoir que dans le sch ma MC 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 MC pr sente encore des oscillations alors que la solution MD para t raisonnable partir de t 0 02 d apr s la figure 2 7 toutes les temp ratures de la solution MC sont positives et le restent bien que la solution MD paraisse un peu plus pr cise Nous suivrons l opinion de Gresho et al 3 en disant que les oscillations du sch ma MC 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 On observera en effet sur la figure 2 8 que la solution MD 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 MC indiquent que ce sch ma ne v rifie pas le principe du maximum que v rifie la solution exacte Supposons par exemple que la variable soit la concentration d une substance dans le domaine d
42. lise l appel des diff rents blocs fonctionnels de MEFISTO est PPMEF OBJ 22 La cr ation du module ex cutable MEF CV2D est r alis e par l instruction LINK EXE MEFCV2D PPMEF LMCV2D INC ELEMLB L MEF L LIBCOM L On trouvera ci dessous l ent te du sous programme ELEMLB de la biblioth que CV2D qui pr cise la fa on de d finir les param tres de l quation l aide des blocs MILL FORC PRND 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 MEFG dans la suite Des probl mes identiques apparaissent avec la m thode des diff rences finies centr es MDFC 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 T U D O lt x lt L U gt 0 dx dx 2 1 avec les conditions aux limites suivantes 1 Oj T et 2 2 T x L T amp T T 2 3 La solution exacte de l quation pr c dente est pie k T x zi T l e EL Poe nu CEEE 2 4 e 7 T Lee
43. llustrer le concept de stabilit nous nous basons sur le syst me d quations diff rentielles ordinaires suivant MT 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 La 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 TRLC et TRLV du programme MEFISTO est le suivant Pere apa REP M OAK AT KT F t 641 2 14 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 sch ma explicite 8 1 sch ma implicite pur 8 1 2 sch ma de Crank Nicholson 6 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 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 X d finis par K A M X 0 i 1 n 2 16 La matrice X est orthonorm e et v rifie X MX 8 l J LJ X K X A 6 pas de sommation 2 17 Dans la base des vecteurs propres le syst me 2 13 s explicite sous la forme des n
44. 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 2 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 ASME AMD 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 ASME AMD 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 T 70 CTL TR 1 8 X L Figure 2 1 Solution de l quation de diffusion convection pour diff rentes valeurs de PE U L D 9 6 p 2 3 4 9 6 7 B 9 l 67 bauten exacte e 6 wt O ay C 1 2 pai a 0 7 6 b 73 o Gu 1 0 o an ey 2 0 _t B 1 B 2 8 3 8 4 8 lambda dt Figure 2 2 Facteur d amplication sch ma d Euler gt Ot any a aa 5 8 OS
45. n le principe du maximum 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 At 42 l ment lin aire et matrice masse consistante si D At h 1 686 _ e D A t element quadratique et matrice masse consistante si 1 400 s lt 1 106 ho Remarque le fait de ne pas v rifier le principe du maximum discret ne signifie pas que la solution oscillera syst matiquement Tout d pend en effet du second membre gui depend 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 a imposer une solution initiale de la forme T H Il 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 maximum 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 Remarque le fait pour un l ment de v rifier le principe du maximum discret ne saurait tre consid r comme une assurance sur la qualit de
46. n pure 0 est repr sent sur la figure 2 5 en fonction de wAt pour diff rentes valeur de 6 On remarque qu un sch ma implicite avec 6 2 gt 0 5 r duit toujours la fr quence semi discr tis e w Comme 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 6 gt 0 5 Remarque l expression 2 42 montre que l int gration explicite O dans le cas de la convection pure D 0 conduit un sch ma inconditionnellement instable On trouve en effet que dans ce cas D DbtoW1 w At gt 1 alors que le sch ma ne peut tre stable que si D gt 0 Une mesure de la dissipation num rique introduite par la discr tisation temporelle est donn e par le rapport t w On 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 u w en fonction de kh s tendent qualitativement au cas bidimensionnel La 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
47. olution exacte de 2 46 aux solutions num riques obtenues en r solvant de fa on exacte le syst me 2 50 dans les cas MD et MC La solution exacte de 2 46 est sinnnx ty e 2 51 Mis 2 T x t 1 x nl i n n La solution exacte du syst me 2 50 est N J 1 sinna Ea t T t 1 a sin Jnae J l N 1 N 1 2 l cosna 2 52 ou a n N 1 et les valeurs propres de K 1M sont MC _ 9 1 cosna A ON EE he pour le sch ma MC 2 53 MD _ 2 2 N41 1 cosna pour le sch ma MD 2 54 Remarque on peut retrouver les valeurs propres du syst me semi discr tis partir du rapport calcul pr c demment cf 2 2 2 On a en effet trouv 2 1 coskh kh 1 2r coskh 1 2 55 SN ERIN I 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 nf n l soit k n dans l expression 2 55 Le param tre d amortissement exact est amp Dk soit n 2 dans le cas pr sent En tenant compte du fait que kh nf N 1 ng et h 1 N 1 on trouve Z 2 N 1 1 cosna 1 2r 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 MC amplifie le param tre d amortissement exact n 92 et que le sch ma MD le r duit dans un rapport sensiblement identique MC MD De fait la va
48. ortantes 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 PECLET NUMERIQUE 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 I 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 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 est loin d tre imp rative dans tous les cas 45 La fo
49. pour suivre cette couche limite Un nombre de P clet num rique P trop lev Des 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 a la bonne vitesse Le rem de consiste alors a 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 comme 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 comme moins imp
50. ra Se SE m EE CURE RE RE re 3 2 2 0 1 8 B e 1 2 2 0 3 2 4 2 5 8 8 2 Figure 2 11 Maillage grossier 55 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 do 0 1 ay 0 5 ay 0 1 isothermes de 0 1 0 9 par pas de 0 1 T odn 0o Figure 2 16 Solution stationnaire maillage grossier do 0 01 ay 0 1 ar 0 02 isothermes de 0 1 0 9 par pas de 0 1 Figure 2 17 Solution stationnaire maillage grossier do 0 01 ay 0 1 ay 0 02 isothermes de 0 1 0 9 par pas de 0 1 2I Figure 2 18 Solution stationnaire maillage fin do 0 01 ay 0 1 ar 0 02 isothermes de 0 1 0 9 par pas de 0 1 Figure 2 19 Solution stationnaire maillage fin do 0 01 ay 0 1 r 0 02 isothermes de 0 1 a 0 9 par pas de 0 1 Figure 2 20 Comparaison des solutions stationnaires obtenues avec le maillage grossier et le maillage fin 58 Figure 2 2 l a Solution transitoire maillage grossier do 0 1 ay 0 5 ay 0 15 t 0 05 isothermes de 0 1 0 9 par pas de 0 1 Figure 2 21 C d t 3 05 r alisation service reprographie du BRGM 86 SGN 194 EA
51. re de dimension de l espace est ng ce qui permet de donner un aspect tout fait g n ral la formulation Soit V un domaine ouvert de Rd et S sa fronti re Soit x xi i 1 2 Ng Un point de V et n n 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 US c q 1 5 b sS S o et Sg 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 est l ensemble vide Nous utiliserons la convention de sommation implicite sur les indices r p t s par exemple si n 2 alors qjnj qin qhn2 Une virgule est utilis e pour indiquer une d rivtion partielle par exemple qi i 0q 0x Consid rons une triangulation du domaine V en l ments VE e 1 2 Ne o ne est le nombre d l ments La fronti re de chaque l ment est not e S et nous supposons V ve 1 6 SES L ensemble U VE repr sente les int rieurs des l ments hors fronti res et l ensemble Sint U S S repr sente les fronti res int rieures de la triangulation L quation 1 1 fait appara tre 3 types de flux le flux convectif g U C 1 7 L l e flux diffusif qe DC 1 8 le flux total g a 1 9 L quation 1 1 peut compte tenu de ces notations s crire sous la forme q f Q waC at 1 10 Les notations suivante
52. rme de l quation de diffusion convection utilis e dans la biblioth que CV2D permet de g n rer un sch ma upwind manuel en contr lant la diffusion artificielle introduite En effet le nombre de P clet num rique dans la direction de l coulement est Vah n a d Fa IV Il suffit donc d augmenter la dispersivit longitudinale pour diminuer P 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 MEFG utilisant une matrice de masse constante est bien meilleure que celle des m mes sch mas utilisant une matrice masse diagonalis e ce qui explique que les r sultats obtenus par la MEFG sont souvent bien meilleurs que ceux obtenus par les diff rences finies centr es pour un sch ma d int gration temporelle du type Euler L utilisation des sch mas MD se justifie lorsque le sch ma d int gration temporelle est explicite Dans la mesure o les blocs TRLC et TRLV du programme MEFISTO ne sont pas con us pour utiliser les avantages d une r solution explicite l emploi de la matrice d
53. 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 En 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 pade FN 2 26 La fr quence w est une fonction r elle du nombre d onde k A est l amplitude de l onde et 7 r La vitesse de phase qui est diff rente de la vitesse d onde en g n ral est 2 27 Cy Il x fF 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 me pour tous les nombres d onde k les modes seront propag s des vitesses diff rentes d o le ph nom ne de dispersion La quantit kx wt est la phase 29 k la longueur d onde et 24 w la p riode Pour un syst me non dispersif la vitesse de phase est une constante Co et w Cok la solution 2 26 est alors celle de l quation monodimensionnelle hyperbolique On peut alors adopter la d finition classique la fonction 2 26 est dispersive si d w dk 0 Dans la plupart des cas l quation traiter est discr tis e par rapport aux variabl
54. rons comme 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 Le champ de vitesse et le param tre physique w 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 ciet num rique et suffisamment faible pour qu il n y ait aucune oscillation dans la solution calcul e Diminuons maintenant en gardant le m me maillage les dispersivit s et le coefficient de diffusion a 0 1 a 0 02 d 0 01 e Le nombre de P clet est maintenant largement sup rieur 2 sur la majeure partie du domaine On 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
55. s causes des oscillations qu on observe parfois dans les solutions num riques Un exemple d utilisation de la biblioth que CV2D 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 o QL OT probleme s crivent aC gt gt gt gt w2 4 Go F DFc e 1 1 concentration porosit cin matique du milieu vitesse de Darcy de l coulement terme source tenseur de diffusion dispersion Ce tenseur est lui m me la somme de deux tenseurs Do do d tant le coefficient de diffusion mol culaire et ll le tenseur identit D 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 En g om trie plane et dans un rep re li l coulement De s crit be a Fo ir O L D a est la dispersivit longitudinale dans le sens de l coulement est la dispersivit transversale dans le plan perpendiculaire l coulement est le module du vecteur vitesse Si U et Uy sont les composantes du vecteur vitesse dans le repere du les composantes du tenseur de diffusion dispersion dans ce rep re 3 2 2 a US GU L T 3 X
56. s seront utiles dans la suite V V n q q n gqg q n qg q n 1 11 n L i n l L r l l ri l l Compte 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 surs 1 12 o g est une fonction donn e d finie sur So La portion So de la fronti re S est une fronti re concentration impos e Il reste d finir une condition sur Sq Nous envisagerons 2 types de conditions qd l sur Se 1 13 q h sur S 1 14 ou 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 Remarque 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 Cap f WdV 0 1 15 V Les fonction d interpolation sont suppos es satisfaire C g surs et les fonctions de pond ration M Deur 1 16 ce qu
57. sse calcul e est alors m wN dV m 0 1 IL ye i ce qui revient affecter au coefficient mj la somme des termes de la ligne i puisque pa N pour les fonctions d interpolation utilis es 5 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 vi f Q N dV Q 2 N E n det CJ En DEP 1 5 8 CALCUL DES SOLLICITATIONS SUR LES FRONTIERES Le flux h DjjC j nj est suppos variable sur la fronti re Les nj 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 amp sur l ar te de r f rence 1 h N E d EV dE l Aa Jo D V ax 08 7 33 08 Nous consid rerons ici le cas d une ar te lin aire ou quadratique ar te de r f rence En i 2 0 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 const
58. sse diagonale et r 1 6 pour le cas MC matrice masse consistente c est dire issue de l application consistente de la formulation de Galerkine On suppose que la solution de 2 30 peut se mettre sous la forme ee 2 31 33 En injectant 2 31 dans 2 30 on obtient U D 1 rL S o t ed ped Ta oa a U SOLES tS S _ LS 2h I J ha y h t rL S a puisque le terme de gauche ne d pend pas de t ce qui donne les deux expressions suivantes B l rL S tr T LS 0 2 32 J h gril j 1 h J b t Bo t 0 2 33 Injectons la condition initiale 2 29 dans 2 32 pour calculer B O A ikx A Lk x h LLRX ikix h oe LS TLe J 2e e J ne Pl ue A k Te 2 cos kh 2 S S LT re Ik h ikh jtl Re lt ae A ikx 21Te J sin kh On en d duit que 8 est donn par l expression iUsinkh h 2D 1 coskh h un 1 2r coskh 1 soit iwsinkhikh 28 1 coskh kh 7 1 2r coskh 1 2 34 ou w ku et Dk sont la fr quence et le param tre d amortissement exacts vus pr c demment La solution exacte de l quation 2 33 est b t e PE 2 35 34 et par analogie avec la solution de l quation initiale 2 28 o O t e itv amp z 2 36 nous r crivons l quation 2 35 sous la forme t iwilet P 2 37 I r sulte de l quation 2 34 que le rapport de la fr quence semi discr tis e w la fr quence exact
59. ti res Sq Remarque En g om trie plane x y dV dx dy 1 dS ds 1 En g om trie axisym trique r z dV 2fr dr dz dS 2fr ds ou s est l abscisse curviligne le long de la fronti re 12 1 5 CALCUL EFFECTIF DES MATRICES ELEMENTAIRES La biblioth que CV2D 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 n 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 N En x 1 31 J N En y o x 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 1 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
60. tte quation montre clairement que pour Ph gt 2 des oscillations qui n existent pas dans la solution r elle 2 4 vont apparaitre L aspect de ces oscillations d pend de la parit de N et leur amplitude de la valeur Pp Il existe de nombreux sch mas de discr tisation appel s upwind permettant de supprimer les oscillations lorsque P est gt 2 Ces 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 a d centrer vers l amont le terme convectif Ainsi le sch ma upwind pur s crit pour u 20 T 2T T 2 8 L L i Pn T T l l 1 ted 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 22 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 MDFC et la MEFG lin aire Pn ee yee et 2 9 Cette quation est identique 2 5 si ce n est que le nombre de P clet de maille effectif P est Pn 2 tanh Pn 2 2 10 Si Pa est remplac par Pp dans l quation 2 5
Download Pdf Manuals
Related Search