Soutenance mémoire
Télécharger le fichier pdf d’un mémoire de fin d’études
Equation de Richards
On cherche a decrire un ecoulement diphasique non miscible eau-air (phases notees respectivement w et a par la suite) dans un milieu poreux. On se place dans un do-maine R2, entierement constitue d’un milieu poreux tel que present ci-dessus. L’equation de conservation de la masse dans (que l’on suppose ferme, sans echange avec l’exterieur), appliquee a chaque phase 2 fw; ag, decrit la dynamique de l’ecoulement : @t( ’s ) + r ( v ) = 0: (1.1)
Les dependances spatiales et temporelles sont ici omises par souci de clarte. Les vitesses v de chaque phase s’expriment a l’aide de la loi de Darcy generalisee [28] : v = K k (s ) [r (p + gz)] : (1.2)
Notons en particulier que, si la direction du gradient de pression n’est pas une des direc-tions principales du tenseur K, alors la direction d’ecoulement du uide et la direction du gradient ne sont pas colineaires.
En l’etat, il s’agit d’un systeme de deux equations a quatre inconnues : pw, pa, sw et sa. Deux contraintes sont necessaires pour fermer ce systeme. La premiere decoule directement de la de nition des saturations : sa + sw = 1:
Nous avons donc une seule inconnue de saturation s := sw, l’autre se deduisant de la formule precedente. La deuxieme contrainte traduit le fait que l’on sait exprimer la di erence des deux pressions inconnues par : pa pw = p:
La fonction p ainsi que ses dependances seront precisees plus loin. On fait alors l’hypothese, classique en hydrogeologie, que l’air circule librement dans le milieu poreux : sa viscosite a est consideree nulle (de l’ordre de 10 6P a:s en realite), si bien que, d’apres la loi de Darcy (1.2) appliquee a l’air, on a : soit ka(s) = 0, ce qui implique que le milieu est entierement sature en eau. Il y a alors degenerescence du systeme en deux equations de Darcy; soit r (pa + agz) = 0. La repartition de la pression d’air est alors hydrostatique :
pa = p0 agz, ou p0 est la pression de l’air a la surface (z = 0), c’est-a-dire la pres-sion atmospherique, que l’on considere constante durant les temps d’observation.
Methodes numeriques pour les ecoulements en milieu poreux
On presente ici rapidement le c ur de ce travail de these. La modelisation des ecoule-ments en milieu poreux necessite l’utilisation d’un schema precis et robuste vis-a-vis de maillages generaux et des non-linearites presentes dans l’equation (1.4). Notre choix se porte sur un schema DDFV (acronyme de Discrete Duality Finite Volume) pour la discretisation en espace, et le schema BDF2 (pour Backward Di erentiation Formula of order 2 ) pour la discretisation en temps. On cherche ensuite a reduire le temps CPU en adaptant le pas de temps et le nombre d’iterations de linearisation au cours de la simulation a l’aide de la theorie des estimations a posteriori.
Motivations
Les contraintes precedentes nous poussent a utiliser un schema de discretisation choisi :conservatif : on travaille avec la forme conservative de l’equation de Richards; precis et robuste dans le sens donne dans la section 1.1. En particulier, il doit ^etre valable sur des maillages generaux sans degradation de l’ordre de convergence, adapte a la discretisation de ux non lineaires et anisotropes, et capable de traiter des heterogeneites, voire des discontinuites spatiales des proprietes physiques.
En n, comme on l’a vu dans la sous-section 1.2.3, les lois constitutives et k sont non lineaires. On doit donc traiter un probleme d’evolution en temps, qui necessite des linearisations pour se ramener a la resolution de systemes lineaires a chaque iteration en temps. Le modele requiert donc des moyens de calcul e caces, une contrainte d’autant plus pertinente si l’on souhaite, par exemple, quanti er les incertitudes de modelisation. Il est particulierement interessant dans ce cas de questionner l’optimalite du pas de temps choisi ainsi que du critere d’arr^et de la procedure de linearisation.
Discretisation du ux di usif
On synthetise ici quelques schemas numeriques utilises pour la discretisation de ux di usifs non lineaires, notamment celui present dans l’equation de Richards.
Parmi les schemas de type elements nis, Knabner [55] a developp une methode d’elements nis mixtes (MFE) pour la discretisation de l’equation de Richards, et parvient a un systeme dont la matrice est de nie positive. Cependant, la mise en place est peu aisee et le surco^ut de calcul par rapport aux methodes d’elements nis classiques, non negligeable.
Les methodes de Galerkin discontinues (DG) forment une classe de schemas popu-laires pour bon nombre de problemes d’applications industrielles. Ce sont techniquement des methodes d’elements nis, mais qui bene cient egalement de proprietes normalement propres aux volumes nis (conservativit locale due a un calcul des ux aux interfaces du maillage, exibilit dans l’utilisation de maillages non conformes). En particulier, une discretisation de l’equation de Richards utilisant la version a penalisation interieure symetrique de ces methodes (SIPG) est consideree dans [73]. Dans la famille des methodes volumes nis, un premier choix «simple» est le schema a deux points [39] (egalement appel TPFA ou VF4), pour lequel on se donne une inconnue locale par volume de contr^ole. Nous verrons cependant dans le chapitre 2 qu’il n’est pas adapte a la modelisation en milieu heterogene et anisotrope.
Le schema diamant [24] s’appuie sur le maillage eponyme, construit en reliant les som-mets d’une ar^ete donnee du maillage de depart aux centres des mailles qui se partagent cette ar^ete. Les valeurs de la fonction inconnue aux sommets de l’ar^ete sont interpolees a partir de ses valeurs aux centres des mailles associees par la methode des moindres carres. Ces quatre valeurs servent ensuite a de nir les ux numeriques. Le schema est valable sur maillages generaux, mais sa construction rend l’analyse delicate. En par-ticulier, les estimations d’erreur de [24] dependent d’une coercivit obtenue sous des conditions geometriques exprimees sur le maillage diamant. De plus, le systeme lineaire associe n’est pas symetrique. Il reste que le schema est regulierement utilise car robuste, et permet de monter en ordre assez facilement.
Les schemas MPFA [1] utilisent des inconnues intermediaires localisees aux centres des ar^etes pour reconstruire un gradient discret sur des sous-volumes du maillage. La solution est supposee a ne sur chacun de ces sous-volumes, et les inconnues intermedi-aires sont eliminees en exprimant la continuite de la solution ainsi que la conservativit des ux. Le choix d’expression de ces equations conduit a di erentes versions : ce sont entre autres les schemas O [1],L [3],G [4] et U [2]. On obtient au nal des ux numeriques consistants, conservatifs et leur expression fait seulement intervenir les inconnues asso-ciees a chaque volume de contr^ole. Le principal inconvenient a retenir est la perte de coercivit et/ou de monotonie dans certains cas (fortes heterogeneites, anisotropie).
Les schemas mimetiques MFD [16] se resument a la recherche d’un produit scalaire discret veri ant certaines proprietes de coercivit et de consistance. En utilisant la formule de Green, et en construisant un operateur de ux approche en dualite avec la divergence discrete naturelle (ces operateurs discrets miment les operateurs conti-nus), on arrive a un systeme lineaire de type point-selle; la monotonie n’est pas assuree. La methode est cependant inconditionnellement coercive sur maillages generaux. Leur generalisation a des problemes non lineaires n’est pas evidente.
Les schemas SUCCES/SUSHI placent des inconnues scalaires au centre de chaque maille et de chaque ar^ete, permettant d’ecrire un gradient discret. Le schema SUC-CES [40] elimine ensuite les inconnues d’ar^etes en les exprimant comme des combinaisons convexes des inconnues de maille. Le nombre d’inconnues est alors restreint, mais des instabilites numeriques peuvent se produire lorsque le tenseur est discontinu. C’est ainsi qu’est apparu le schema SUSHI [41], qui fait le choix de garder les inconnues d’ar^etes a travers lesquelles le tenseur n’est pas regulier, et de fermer le systeme en ecrivant la continuite des ux discrets a travers ces ar^etes. La precision du schema s’en trouve amelioree, au prix d’un nombre d’inconnues plus elev . Notons qu’une approche uni ee des methodes volumes nis mixtes, MFD et SUCCES a et developpee dans [35]. Les schemas VAG [42] utilisent des inconnues supplementaires situees aux n uds du maillage pour de nir un gradient discret sur des sous-volumes de chaque maille, qui est ensuite stabilise a n d’eviter des degenerescences indesirables. On associe a chaque n ud un volume dependant des proprietes physiques autour de ce n ud, puis on ecrit
une equation exprimant un bilan de ux nul sur ce volume. A ce prix, on obtient un schema inconditionnellement coercif sur maillages generaux. Le schema a et im-plement pour des ecoulements diphasiques incluant l’equation de Richards dans [43]. Notons cependant que le choix de la redistribution des volumes de n ud reste relative-ment empirique.
Les schemas DDFV ont et introduits dans [52]. Leur principe est, comme pour les schemas VAG, d’ajouter des inconnues aux sommets du maillage. Un maillage secondaire est construit a partir du maillage de depart (dit primaire), chacun de ses volumes etant associe a un n ud unique. L’equation continue est alors integree sur chaque volume de ces deux maillages, ce qui conduit apres application de la formule de Green a ap-procher les composantes normale et tangentielle du ux continu. Ces schemas ont et utilises sur des problemes varies ces dernieres annees, allant de la di usion scalaire [34] a l’electrocardiologie [23], entre autres. La symetrie du probleme continu est preservee (sauf dans le cas de conditions de bord mixtes Dirichlet/Neumann, que l’on etudiera au chapitre 2), et les resultats numeriques exhibent une approximation particulierement precise du gradient, comme cela a et note dans [51]. Pour ces raisons, auxquelles s’ajoute sa simplicite d’apprehension et de mise en oeuvre sur des domaines a deux di-mensions, en tant qu’extension «naturelle» du schema a deux points, c’est ce choix que l’on retiendra dans ce manuscrit.
Discretisation du terme instationnaire
L’equation de Richards (1.4) est instationnaire, il est necessaire de discretiser la derivee en temps @t ( ). Ayant choisi un schema d’ordre 2 pour la discretisation du ux di usif, il semble naturel d’avoir la m^eme exigence pour le terme en temps. Remarquons pour commencer qu’il ne nous est pas permis d’utiliser un schema explicite ici, comme cela a et mentionn dans [44]. En e et, la teneur en eau est constante a saturation, donc non inversible, et le schema associe est mal pose. Le schema de Crank-Nicolson est un choix populaire, particulierement adapte a l’analyse a posteriori presentee au chapitre 3 [33]. Il peut neanmoins se reveler instable et introduire des oscillations parasites sur la solution lorsque la condition initiale n’est pas reguliere [65, 70], comme dans certains cas tests presentes dans les chapitres 2 et 4. Il est egalement possible d’utiliser un schema de Runge-Kutta, mais la presence de plusieurs etapes de calcul ralentit l’execution du code, tout particulierement lorsqu’il est couple avec un solveur non lineaire. Il a egalement et relev dans [73] que le schema de Runge-Kutta diagonalement implicite d’ordre 3 (DIRK3) est inutilisable lorsqu’une partie du domaine est saturee. En e et, la formule de nissant la seconde etape du schema est de nouveau mal posee a cause de l’absence de terme implicite. La simulation de cas tests raides (voir notamment le cas test de Polmann, sous-section 2.4.2) nous demande en n d’^etre exigeants quant a la stabilite du schema retenu. Le comportement des schemas en temps sur des problemes raides est bien decrit par l’etude de la A-stabilite [58]. Parmi les methodes lineaires multipas, on rejette donc les methodes d’Adams-Moulton et d’Adams-Bashforth au pro t de la methode BDF2 [25], la seule a ^etre A-stable. Il est par ailleurs notable que la montee en ordre degrade cette A-stabilite : on peut montrer qu’un schema implicite lineaire, multipas et A-stable est au plus d’ordre 2 (il s’agit de la fameuse seconde barriere de Dahlquist, voir [27]).
Estimations a posteriori
La seconde partie de ce travail concerne les estimations a posteriori. Celles-ci visent de maniere generale a obtenir une borne de l’erreur entre la solution exacte du probleme continu et la solution approchee issue de la simulation numerique, et ce a partir de la solution approchee seulement. Le resultat est fort : m^eme si on ne conna^t pas la so-lution exacte (ce qui arrive dans l’immense majorite des cas), on est capable de donner un majorant, voire dans certains cas un minorant de l’erreur commise par la resolu-tion numerique. Lorsque cette borne est entierement calculable, on parle d’estimations garanties. Notre objectif est d’avoir plus precisement un contr^ole sur la distribution de l’erreur, et d’utiliser cette information pour elaborer une strategie d’adaptation du pas de temps et du nombre d’iterations de linearisation, le maillage etant xe par ailleurs; nous avons donc besoin de bornes locales en temps. Le principe general est de choisir a chaque iteration en temps, un pas de temps «optimal», au sens ou il equilibre les sources d’erreur provenant de la discretisation en temps et de la discretisation en espace. On determine aussi un critere d’arr^et pour la procedure de linearisation qui contraint les erreurs de linearisation a ^etre negligeables devant les erreurs en temps et en espace.
Plusieurs strategies sont a notre disposition pour etablir des estimations a posteriori. La theorie a debut avec les travaux de Babuska et Rheinboldt [9] dans les annees 1970, qui ont propose des estimations par residu etendues ensuite par Verfurth [77]. Cette meth-ode borne l’erreur entre les solutions exacte et approchee dans une norme d’energie. Elle utilise une evaluation de residus locaux relatifs a la forme forte de l’equation a resoudre, mais font generalement intervenir une constante multiplicative di cile a im-plementer [79]. L’estimateur de Zienkiewicz-Zhu (voir [80]), donne une estimation locale du gradient de la solution exacte a l’aide de la solution approchee. Cette methode est populaire gr^ace a sa simplicite de mise en oeuvre et son faible co^ut de calcul, mais ne per-met pas de bene cier d’une estimation garantie. Les methodes hierarchiques enrichissent l’espace d’approximation (par ra nement en temps et en espace) pour calculer les esti-mateurs. Elles sont souvent utilisees pour adapter le maillage. La methode par residus equilibres resout des problemes aux limites locaux, ou les donnees de bord sont choisies de maniere a equilibrer le residu interieur. L’article de Bank et Weiser [10] contient des idees fondamentales (notamment l’hypothese de saturation ainsi que l’equilibrage des donnees de bord) reprises par beaucoup de travaux sur les methodes hierarchiques et par residus equilibres. L’ensemble de ces techniques est synthetis dans [6] pour les elements nis.
Concernant les methodes de volumes nis, Ohlberger obtient une estimation de l’erreur en norme L1 pour des equations de convection-di usion non lineaires, respectivement pour des volumes nis centres par maille dans [62] et centres par sommet dans [63]. Une strategie d’adaptation de maillage y est egalement proposee. Plus speci quement pour un schema DDFV, Omnes [64] propose une estimation de l’erreur en norme d’energie, en passant par une formulation variationnelle du schema qui permet d’utiliser des tech-niques bien connues pour les elements nis. Des estimations sont obtenues a la fois sur le maillage primaire et sur le maillage secondaire, et font appel a une decomposition de Helmholtz-Hodge.
Notre travail rentre dans la categorie des methodes de ux equilibres [68], fondees sur la reconstruction de ux continus a partir des ux discrets du schema, equilibres en un certain sens. Ce type de methode a recu une attention particuliere dans de nom-breuses etudes ces dernieres annees : des problemes d’elasticit traites par elements nis ainsi que l’equation de Poisson dans [30, 57], des schemas DG pour une equation de convection-reaction-di usion dans [37], des volumes nis pour les ecoulements mul-tiphasiques dans [32]. Plus recemment, des developpements theoriques ont uni e une classe generale de discretisations en espace, d’abord pour l’equation de la chaleur [38], puis pour un probleme parabolique non lineaire [33]. Dans ce dernier article, des es-timations (bornes superieure et inferieure) robustes et completement calculables sont ecrites a l’aide d’une norme duale espace-temps. D’autres estimations pour l’equation de Richards peuvent ^etre trouvees dans [13], mais utilisent la transformee de Kirchho . On s’interesse pour notre part a la formulation de l’equation de Richards relative a la charge hydraulique, sans application de la transformee de Kirchho ; un schema con-servatif tel que DDFV est justement concu pour travailler avec les variables physiques. Notre travail s’organise en trois parties. On commence par ecrire une borne superieure dans l’esprit de [33, 68]; puis, on propose des reconstructions espace-temps adaptees a la discretisation DDFV-BDF2 de l’equation de Richards, qui necessitent une reformulation du schema BDF2. La non-linearit du terme en temps necessite un traitement parti-culier : on equilibre alors le ux temporel aussi bien que le ux spatial, en coherence avec la norme espace-temps utilisee dans les estimateurs. En n, on distingue plusieurs composantes d’erreur dans notre estimation : une erreur provenant de la discretisation en temps, une erreur provenant de la discretisation en espace et une erreur de linearisa-tion. On propose un algorithme adaptatif qui permet de choisir un critere d’arr^et pour la procedure de linearisation, ainsi qu’un pas de temps qui equilibre les erreurs en temps et en espace. Soulignons que les resultats presentes restent generaux et peuvent s’adapter a d’autres discretisations.
La matrice du systeme (2.17)-(2.18) garde la m^eme structure tout au long de la sim-ulation, et elle doit ^etre assemblee un grand nombre de fois, notamment a cause de la procedure de linearisation. C’est pourquoi on utilise une methode de renumerotation, qui est appliquee une seule fois sur les centres et les sommets des mailles primaires (sauf les sommets du bord de type Dirichlet, dont les mailles secondaires associees ne par-ticipent pas au systeme), en pretraitement. On choisit ici la permutation inversee de Cuthill et McKee [26]. On perd la structure par bloc de la matrice initiale, mais les coe cients non nuls de la nouvelle matrice sont resserres autour de la diagonale (voir la gure 2.8), conduisant a une largeur de bande minimale et a un gain en temps CPU signi catif.
Initialisation de la procedure de linearisation
Ordre 0 Il s’agit de l’initialisation la plus simple : n;0 = n 1.
Ordre 1 On choisit n;0 = P1(tn), ou P1 est la fonction a ne qui interpole en tn 2 et en tn 1, rencontree dans la sous-section 2.2.2. On trouve n;0 = 2 n 1 n 2.
Ordre 2 On choisit n;0 = P2(tn), ce qui donne : n;0 = 3 n 1 3 n 2 + n 3.
Des initialisations d’ordre quelconque d 0 peuvent ainsi ^etre implementees, en utilisant le polyn^ome d’interpolation de Lagrange qui interpole les vecteurs n 1; ; n d 1 en tn 1; ; tn d 1. Les coe cients du polyn^ome servant a l’initialisation sont ceux du triangle de Pascal avec signes alternes : pour une initialisation d’ordre d, on a : n;0 = Pd+1( 1)i+1 d+1 n i.
|
Table des matières
Table des mati`eres
Liste des Figures
Liste des Tableaux
1 Introduction
1.1 Motivations
1.2 Mod`ele physique
1.2.1 Milieu poreux
1.2.2 Equation de Richards ´
1.2.3 Lois de fermeture
1.2.4 Mod´elisation du sous-sol
1.3 M´ethodes num´eriques pour les ´ecoulements en milieu poreux
1.3.1 Motivations
1.3.2 Discr´etisation du flux diffusif
1.3.3 Discr´etisation du terme instationnaire
1.3.4 Estimations a posteriori
1.4 Contributions de la th`ese
1.5 Plan du manuscrit
2 Sch´ema DDFV pour l’´equation de Richards
2.1 Diffusion en milieu h´et´erog`ene anisotrope
2.1.1 Limitations du sch´ema VF4 pour l’´equation de Poisson
2.1.2 Sch´ema DDFV pour la diffusion lin´eaire
2.1.3 Extension `a la diffusion non lin´eaire et aux conditions aux limites mixtes Dirichlet/Neumann
2.2 Equation de Richards ´
2.2.1 Discr´etisation en espace
2.2.2 Discr´etisation en temps
2.2.3 Lin´earisation
2.2.4 Evaluation du tenseur ´
2.3 R´esolution num´erique
2.3.1 Assemblage du syst`eme lin´eaire
2.3.2 R´esolution du syst`eme lin´eaire
2.3.3 Initialisation de la proc´edure de lin´earisation
2.4 Cas tests
2.4.1 Infiltration en milieu homog`ene isotrope : validation (TC1)
2.4.2 Infiltration en milieu homog`ene isotrope : cas raide (TC2)
2.4.3 Quarter five spot en milieu h´et´erog`ene anisotrope (TC3)
3 Estimations a posteriori pour l’´equation de Richards 53
3.1 Probl`eme continu
3.1.1 Hypoth`eses sur les donn´ees
3.1.2 Solution faible, d´efinition du r´esidu
3.2 Sch´emas en temps et lin´earisations
3.2.1 Cadre discret g´en´eral
3.2.2 Sch´emas d’Euler implicite et de Crank-Nicolson
3.2.3 Sch´ema BDF2
3.2.4 Synth`ese
3.3 Estimations a posteriori par flux ´equilibr´es
3.3.1 Notion d’´equilibrage de flux
3.3.2 Lien avec les sch´emas
3.3.3 Estimations pour un solveur non lin´eaire exact
3.3.4 Estimation pour le probl`eme lin´earis´e
3.4 Choix de reconstructions pour le sch´ema DDFV-BDF2
3.4.1 Reconstructions de la charge
3.4.2 Reconstructions de la teneur en eau
3.4.3 Reconstructions du flux spatial
3.4.4 Reconstructions du terme source
4 Algorithme et r´esultats num´eriques
4.1 Algorithme propos´e
4.2 Cas tests
4.2.1 Cas test analytique
4.2.1.1 Espace de reconstruction du flux Hdiv
4.2.1.2 Simulation sans raideur
4.2.1.3 Simulation avec raideur en temps
4.2.2 Cas test de Polmann
4.2.3 Sol h´et´erog`ene `a fronti`ere curviligne
Conclusion
Programmation efficace en MATLAB
Bibliographie
Télécharger le rapport complet
