Détermination et étude numérique de conditions aux limites transparentes. Application à la simulation du procédé CHFR de logging de la société SCHLUMBERGER.

Copyright (c) Kout HASSeNE

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation;

with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts. A copy of the license is included in the section entitled "GNU

Free Documentation License".

Détermination et étude numérique de conditions aux limites

transparentes. Application à la simulation du procédé CHFR de

logging de la société SCHL UMBER CER

.

Hassene KOUT Encadrement :
DESS Ingenierie Mathematique option

calcul scientifique et modélisation me-- A Schlumberger RPC :

canique à l'Université de Savoie. Martin G.Luling,

année universitaire (2002-2003) Isabelle Dubourg.

A l'INRIA de Rocquencourt : Jacques Henry.

Résumé

Dans le domaine de la prospection pétrolière, les mesures des propriétés physiques de la roche doivent être en phase de production c'est le reservoir monitoring.

Malgré les conditions opératoires difficiles ( notamment la présence d'un casing), Schlumberger a développé un outil efficace pour mesurer la résistivité en trou fermé le CHFR, un outil d'injection de courant.

Afin de maItriser et étendre son utilisation, une simulation fiable est souhaitable. Deux méthodes sont actuellement utilisées. Le modéle des télégraphistes approche le problème par une résolution analytique simpliste. Une autre résolution est implémentée dans le programme Cwnlat. Elle résout par méthode d'éléments finis les équations aux dérivées partielles qui gouvernent la distribution du courant.

Au cours d'un log simulé Cwnlat doit résoudre plusieurs fois le même problème, ce qui est lourd en temps de calcul.

L'étude menée dans ce stage a permis d'optimiser le temps de calcul de Cwnlat en appliquant une méthode de factorisation.

La méthode de factorisation permet d'effectuer les calculs sur un sous-domaine fixe qui contient l'ensemble des positions de la sonde en déterminant des conditions aux limites sur le bord du sous-domaine qui résument exactement le comportement de la solution sur le domaine enlevé. Le nombre des inconnues est ainsi réduit.

Une expérimentation sur une configuration simplifiée a montré une bonne optimisation du temps de calcul.

Le rapport fournit des figures qui illustrent l'étude. Le rapport a été rédigé en juillet.

Table des matières

1 Introduction 7

2 Le cadre physique 9

2.1 La vie d'un puits de pétrole 9

2.2 Le modèle pour l'interprétation des mesures 10

2.3 Contexte et but du stage 11

2.4 Modéle adopté dans la simulation 13

2.5 Fonctionnement de l'outil 13

3 Le cadre mathématique 15

3.1 Mise en équations 15

3.1.1 Les conditions aux limites 17

3.2 Méthode des éléments finis 22

3.2.1 Formulation variationnelle 22

3.2.2 Forme discrete 23

3.2.3 Discrétisation Q1 24

4 Méthode de Factorisation 29

4.1 Motivation 29

4.2 Calcul formel de la méthode de factorisation 31

4.3 Forme discréte 34
4.4 Prise en compte des conditions aux limites transparentes dans

Cwnlat-bis 42

4.4.1 Formulation variationnelle 43

4.4.2 Forme discrete 44

4.4.3 Discrétisation Q1 46

4.5 Tests 48

5 Conclusion 59

Bibliographie 60

TABLE DES MATIERES

Chapitre 1

Introduction

Schlumberger Ribout Product Center (SRPC) est un centre d'études et production de la section Wireline & Testing de la branche Oilfield Services de la grande multinationale Schlumberger. Le centre est situé à Clamart au sud de Paris, et emploie environ 500 personnes.

Les activités de SRPC se concentrent autour de la conception et la simulation des outils de mesures dans les puits de pétrole. Le centre dispose en outre de deux puits d'essais pour tester les outils. La chaIne des opérations menées est assez complete on va de la modélisation par CAO à la réalisation, en passant par tous les aspects électroniques, mécaniques, matériaux, etc. Cependant, SRPC a de plus en plus tendance à sous traiter quelques étape de la production.

SRPC est en liaison étroite avec le centre de Schlumberger situé à Houston au Texas, de nombreux cadres se déplacent fréquemment de Paris à Houston et vice versa. Et comme à Houston, le niveau du prix du pétrole influe considérablement sur l'activité du centre.

Parmi les unités de SRPC, le projet CHFR travaille sur l'outil CHFR de mesure de résistivité en puits fermé. Il est constitué d'une dizaine de personnes qui travaillent entre autre sur l'éléctronique et la simulation.

A cause d'un problCme de temps d'éxécution de CWNLAT qui est long, le projet a cherché une collaboration avec le projet ONDES à l'INRIA Rocquencourt, C'est ainsi que s'est présenté le stage.

Le projet ONDES à l'INRIA travaile sur le modélisation et résolution de problCmes de propagation d'ondes électromagnétiques, acoustiques etc. Il comprend environ six personnes.

Le stage était donc dans un contrat de collaboration entre Schlumberger et l'INRIA, il s'est déroulé en partie à Clamart mais surtout à Rocquencourt, et ce de Avril à juillet 2003.

Introduction

Chapitre 2

Le cadre physique

2.1 La vie d'un puits de pétrole

Après avoir localisé (par mesures sismiques) une formation rocheuse susceptible de contenir du pétrole (souvent en-dessous d'une couche de roches imperméables tel l'argile), on procède au forage du puits. Le forage creuse jusqu'un peu en-dessous de la formation fertile, et pendant cette opératoin, une boue spéciale est déversée dans le puits pour retenir le pétrole.

Aprés ce forage, le puits est en"trou ouvert" (open hole). On en profite pour faire une série de mesures des paramètres de la roche. Ces mesures sont réalisées par des instruments électriques et/ou nucléaires, elles vont servir à identifier un modéle physique de la roche et à inspecter sa contenance effective en hydrocarbures.

Une fois que ces mesures ont montré la rentabilité du puits, on vient assurer sa tenue mécanique par un cylindre en acier le casing, qui va jusqu'un peu en-dessous de la formation fertile. Il est installé par morceaux et fait couramment 3000 m de longueur. Le casing est fixé par une cimentation autour de lui. Enfin on procède à des perforations par explosifs au niveau de la formation pour laisser sortir le pétrole.

A partir de ce moment, le puits peut être mis en production, les dispositifs nécessaires sont alors mis à sa surface (well head). Le puits est "fermé", c'est pourquoi on parle de puits en "trou fermé" (cased hole).

Le suivi régulier de l'état du puits reste cependant un élément important de l'exploitation. Des mesures doivent régulièrement être faites pour contrOler la qualité du pétrole qui est produit et la structure de la roche (surtout contenance en pétrole et état mécanique). Le but est de savoir s'il faut arrêter la production dans un puits ou si l'on peut continuer à produire

et pour combien de temps. Cet ensemble d'opérations de prospection pendant la production est connu sous le nom de"reservoir monitoring".

Les phénomènes et grandeurs physiques qui peuvent déterminer la vie d'un puits sont nombreux. De l'infiltration d'eau dans la roche jusqu'à l'influence d'autre puits voisins, des réponses précises et fiables aux questions posées sont nécessaires pour faire un bon reservoir monitoring.

Les contraintes de coüts et de possibilités techniques obligent à"faire parler la roche " a travers le casing toutes les mesures envisageables sont filtrées par cette géante électrode. Ceci rend les mesures surtout électriques, trés difficiles.

De plus, il faut arrêter la production pendant la mesure pour pouvoir opérer à l'intérieur du casing, une boue spéciale remplit ce dernier pour contrecarrer la pression du pétrole. Ceci implique que les mesures doivent être faites le plus rapidement possible. D'autre part, l'usure rapide des instruments de mesures qui descendent dans le puits est rapide, il est difficile donc de recommencer la mesure en cas d'erreur.

En résumé, un bon reservoir monitoring nécessite des outils de mesures fiables et rapides.

Actuellement, pour la mesure de résistivité de la roche, deux types de mesures sont faites des mesures nucléaires où l'on envoie des radiations et recueille les échos, des mesures par injection de courant en profondeur à travers la casing et on le recueille au well head.

En règle générale, ces trois procédés fournissent ce que l'on appelle un log le tracé d'un signal de sortie en fonction de la profondeur dans la zone fertile. En fait l'outil de mesure fonctionne par station la mesure est répétée à plusieurs profondeurs du puits.

L'étape suivante est de faire la liaison entre la resultat de la mesure et les paramètres physiques intéressant de la roche (tel sa densité, sa résistivité, sa nature, etc.) qui permettront de tirer des conclusions en matière de réservoire monitoring. C'est le métier d'interprétation.

2.2 Le modèle pour l'interprétation des mesures

L'interprétation des mesures n'a pas de sens si l'on ne se donne pas un modèle de la roche, ce qui revient en fait à faire un choix des zones dans les quelles on suppose les paramètres physiques constants.

La première hypothèse serait de dire que la roche est homogène par couches, ainsi un paramètre physique sera attribué à une profondeur et on a

un modèle unidimensionnel à symétrie de révolution autour du casing. Cependant, dans une même couche, les paramètres physiques peuvent varier radialement par rapport au casing. c'est le cas par exemple des zones limitrophes au casing qui on été influencées par le contact avec la boue de forage ou avec le ciment de fixation, on parle couramment de zones envahies. Si l'on veut rendre compte de ces effets, il faut considérer un modèle bidimensionnel (dépendance de la profondeur et du rayon), en première approximation avec touj ours une symétrie de révolution autour du casing. Le modèle complet sera tridimensionnel, où l'on ne suppose plus la symétrie de révolution. ce dernier modéle est indispensable pour rendre compte par exemple des infiltrations d'eau d'une seule direction.

Le modèle est élaboré en phase de trou ouvert, où les mesures précises et fiables permettent de comprendre la roche. Cependant, les mesures ne peuvent être faites qu'avant l'installation du casing et donc ne peuvent être répétées pour suivre l'évolution de la roche. Actuellement, le dernier outil d'injection en trou ouvert ( le HRLA ) permet, en ayant plusieurs profondeurs d'investigation de constituer un modèle bidimensionnel de la roche en matière de résistivité.

Une fois constitué, le modèle en trou ouvert sera une référence pour les mesures ultérieures. Le reservoir monitoring consiste en fait à déceler les chagements d'état de la roche par rapport à ce qui était avant l'exploitation mais en adoptant toujours le même modèle.

2.3 Contexte et but du stage

La simulation du procédé de mesure est un outil trés apprécié en reservoir monitoring. Elle permet de prévoir plus ou moins la réponse de l'outil et permet de décider de son utilité dans une situation précise. Si l'on veut localiser ou mesurer un phénomène dans la roche ( par exemple une infiltration d'eau ), il serait profitable de vérifier au préalable par simulation si l'outil sera sensible à ce phénomène, et si oui dans quelles mesures.

La simulation d'un procédé de mesure nécessite en premier lieu une modélisation de la roche, ce qui a été fait dans les mesures en trou ouvert. Ensuite, il faut modéliser et mettre en équation intégrales et/ou différentielles les phénomènes physiques qui entrent en jeu. Il faut enfin un outil informatique pour résoudre les équations et calculer les grandeurs qu'aurait sorties l'outil en réalité. C'est en examinant le résultat des simulations que

l'on va décider de l'utilisation ou non d'un outil particulier pour mesurer ou quantifier un certain aspect des propriétés physiques de la roche.

Ce stage s'inscrit dans les efforts menés à Schlumberger pour la simulation du CHFR ( Cased Hole Formation Resistivity ) qui est un outil d'injection de courant et qui sert à mesurer la résistivité dans différentes parties de la roche en trou fermé.

La mesure de résistivité est l'une des plus anciennes expertises de Schlumberger dans les service pétroliers. La résistivité est une information précieuse dans la prospection pétrolière puisqu'elle indique directement la saturation de la roche en hydrocarbures. Elle permet donc au départ de voir la rentabilité du puits et pendant le reservoir monitoring de suivre à quelle vitesse la roche perd son pétrole, et donc d'estimer l'espérance de vie du puits.

Le CHFR est une sonde qui descend au fond du puits et puis remonte en s'arrêtant à des stations pour injecter latéralement un courant à travers le casing. Le courant stabilisé à une valeur fixe retourne à la surface du casing, tandis que des mesures de potentiel en-dessous du point d'injection permettent d'estimer le courant qui pénètre dans la formation. Le courant et le potentiel permettent d'estimer le résistivité de la formation.

Le CHFR est actuellement commercialisé pour le reservoir monitiring et il a été utilisé plusieurs fois à travers le monde, cependant on ne dispose pas d'une simulation fiable de son comportement. Une bonne simulation sera surtout utile pour adapter l'outil à de nouvelles situations, tel par exemple une avance d'un front d'eau sur le puits ou une élévation du niveau de la nappe phréatique en dessous.

Schlumberger dispose d'un code de calcul Cwnlat qui a été adapté au cas du CHFR, et d'une résolution analytique du phénomène électromagnétique dans la roche.

Le but de ce stage est d'essayer d'optimiser le temps de calcul de ce code de calcul Cwnlat, en appliquant une méthode de factorisation de domaine qui sera présenter et tester plus loin.

2.4 Modéle adopté dans la simulation

Le modéle adopté pour la simulation de CHFR dans ce stage sera bidimensionnel, avec symétrie de révolution autour du casing. Il est résumé dans la figure ci dessous, qui est une coupe du domaine complet et qui suffit à la simulation puisqu'on a supposé une symétrie de révolution autour du casing. On suppose de plus dans le modéle que la formation s'étend à l'infini, ce qui n'est pas une trés grave approximation.

2.5 Fonctionnement de l'outil

Le courant est injecté à la surface intérieure du casing, et il retourne à la surface comme le montre le figure 2. du fait de le forte conductivité de l'acier du casing, la plus grande partie du courant cirdule dedans et seule une faible partie traverse la formation avant de regagner le Well head. Les lignes de courant ( tangentes au vecteur densité de courant ) sont donc surtout concentrées dans le casing, Elles partent toutes du point d'injection et reviennent au Well head. On note J le courant total dans le casing et j le courant qui fuit normallement par unité de profondeur à la surface extérieure de casing. Le contraste de résistivité entre le casing et le reste du domaine (environ 10⁷) est tellement grand que le casing est une équipotentielle. De ce faits, les lignes de courant à sa sortie sont quasi normales et on a j=-dJ/dz. Pour plus de détaille voir la référence.

Environ cinq métres en dessous de l'injection, trois potentiométres mesurent le potentiel électrique V aux points A,B et C. L'outil calcule alors la dérivée seconde du potentiel par rapport à la profondeur qui, divisée par Rc, n'est autre que j au point B. En pratique la mesure est un peu plus compliquée vu qu'il faut estimer Rc par une premiére phase de calibrage.

L'outil fait cette procédure à plusieurs stations en remontant tout au long de la zone qui nous intéresse de la formation ( de 300 à 800 m d'epaisseur ). L'hypothése d'interprétation consiste à dire que la résistivité de la roche en face du point de mesure B (et non celle du ciment ou de la zone envahie) est proportionnelle au quotient VB/jB. Le résultat final des mesures est donc k.VB/jB le long de la formation, le coefficient k est une constante estimée d'après la géométrie du casing.

Pour faciliter le traitement électronique, l'outil injecte un courant

alternatif (AC). Cependant, le procédé est limité à une trés basse fréquence (quelques Hertz) de peur que l'épaisseur de peau dans le casing ne soit trop petit et fasse un écran total entre le courant injecté et le formation.

^-*E)=-,i.^?

-*

= -i.W.,i.H (3.5)

-* rot(

-* H ?t

Chapitre 3

Le cadre mathématique

3.1 Mise en equations

champ électrique ^-E, le champ magnétique

Le phénomène physique de circulation de courant qui se manifeste
dans le fonctionnement de l'outil obéit aux équations de Maxwell. Avec le
^**-H, la permittivité ,i, la conduc-

tivité électrique or, la permittivité électrique å0 et le densité volumique de courant ^*-I, ces équations sont :

-* rot(

^-*E) = -,i. ?-* ?t (3.1)

H

-* rot(

-* H) =

*-I +E0. ?-* ?t (3.2)

E

div(^-* E ) = 0 (3.4)

Le courant injecté étant sinusoIdal à la fréquence f = _2ð, on se

W

place en régime harmonique où l'opérateur de dérivation se remplace par la multiplication par i.W (on adopte cette convention et non -i.W, mais c'est équivalent) avec i² = -1. Seules les trois premières équations sont utiles à la modélisation et elles donnent :

L'équation(3.6) nous donne :

Dans le modèle qu'on a adopté pour la simulation au chapitre précédent, on suppose une symétrie de révolution autour de l'axe z (l'axe du casing) en ce qui concerne les propriétés physiques des différents milieux du domaine (c'est à dire or et ii). Il est naturel alors de supposer une telle symétrie dans les champs électriques et magnétiques.

Dans ce cas (mode transverse magnétique - TM) ces champs s'expriment par :

*-- H = H(r, z).^--* u9(3.8)

* --E = Er(r, z).^--* ur + Ez.^--* uz (3.9)

Q uelques valeurs numériques :

Dans la roche or est comprise entre 0.1 ^Ù-1m-1 et 100 Ù^-1m^-1, ii = 4.ð.10^-7

5.I

Dans le casing or 5.10⁶Ù^-1m^-1, ii = 100.4.ð.10^-7 5.I

Dans le domaine C0 = 36.ð.1095.I, et enfin f vaut quelques Hertz.

1

On constate déjà que le produit ù.C0 est complètement négligeable devant or, ce qui permet de le supprimer. L'équation devient alors avec la seule inconnue H, la composante orthoradiale du champ magnétique, dépendant de r et de z :

L'équation (3.11) est une équation aux dérivées partielles en H, les propriétés sur le courant lui fourniront les conditions aux limites nécessaires

pour avoir un problème bien posé dans un certain domaine géométrique I. La simulation du CHFR se résume donc à la résolution de ce problème aux dérivées partielles avec méthodes numériques adéquates.

3.1.1 Les conditions aux limites

On résout donc l'équation (3.11) dans un domaine tridimensionnel I de frontière 8I. Dans le cas du CHFR, I est le demi espace en dessous de la surface, privé de l'intérieur du casing (la boue qui le remplit est isolante). On choisit aussi d'enlever toute la partie dans la continuation du casing, ceci correspond en fait à une simplification qu'on verra plus loin. Une coupe du domaine par symétrie de révolution autour de l'axe, est donnée dans la figure 3.1, comme adopté dans le chapitre précédent, ce domaine s'étend à +8 en r et à -8 en z. Le bord 8I est donc la réunion des différents segments Fi,

Fr, F1, F0, F2 et F.

Pour mettre en évidence les propriétés physiques qui conduisent aux conditions aux limites pour notre équation (3.11), on utilise le théorème d'Ampère, qui dit que la circulation du champ magnétique le long d'une courbe fermée vaut le flux total du courant à travers toute surface limitée par cette courbe.

On prend un disque horizontal D centré sur l'axe du casing et de rayon r1 et dont le centre est au-dessus du casing shoe. Vu que la boue est isolante, le seul courant circulant à l'intérieur de D est éventuellement celui du cable. Au dessus de l'électrode d'injection, ce courant vaut -I0 car le sens conventionnel est celui de l'axe des z. En dessous de l'électrode d'injection, ce courant est nul. Une première approximation serait de dire que le courant est nul aussi en dessous du casing shoe bien que, la roche étant intacte dans cette zone, on puisse avoir des courants. Mais cette approximation revient à dire que le courant qui y circule est très faible car de façon prépondérante il remonte vers la surface pour rejoindre l'électrode de recueil.

?-

Ainsi, commeHH.^-? n= 2.ð.r1.H(r1,z), on a comme conditions

?D

aux limites sur sur F1 :

2.ð.r1.H(r1,z) = -I0 (3.12)

Le cadre mathematique

z

1'

r

1'

Ù

r

1'

1

r

2

r

1

1'

i

L

1'

1'

0

1'

s

1'

2

1'

Axe de symetrie de revolution

FIG. 3.1 Coupe du domaine de resolution, le casing hachuré. Les points noirs délimitant les frontières designees.

sur F0 U F2 on a la condition aux limites suivante :

2.ð.r1.H(r1, z) = 0 (3.13)

En ce qui concerne la frontière F, où z = z_ ou r = r+ , on a un champ magnétique nul car très éloigné des sources. Ainsi on aboutit aux conditions aux limites :

sur F : 2.ð.r.H(r, z surface) = 0 (3.14)

Restent les électrodes d'injection et de retour, qui sont respectivement un anneau cylindrique petit en épaisseur (1 cm environ) et l'électrode de retour = anneau à la surface. Sur ces zones, le plus simple est de sup-poser un potentiel électrique constant, ou un champ électrique normal à la surface de l'électrode, c'est équivalent. En notant ^*--m la normale extérieure et

*-- r(.)
· ^*--m la dérivée normale, on a donc sur Fi *-- m= ----* uret

8 8m

On va supposer sur les électrodes d'injection et de retour des conditions aux limites plus générale du type :

1

sur Fi et Fr : .(--* m A * ^--u9)
· ^--*rot(H.^--* u9) -- 'y.H = C (3.16)
a

Ici G est un champ électrique imposé et 'y est en quelque sorte l'impédance de l'électrode. Cette condition aux limites est valable sur la partie bidimensionnelle de 8I correspondante.

Finalement, en utilisant l'expression du rotationnel en cooordonnées

--* 8H

cylindriques r et z :rot(H.--*uz ₊18r.H

u9) = -- 8z .--* 8r .--*

.ur et en posant u=r.H,

r

on obtient de l'équation (3.11) :

8r( ¹

8 a.r^.8u 8r ) + 8 8z ( 1a.r.8u8z ) -- i.ù.u _r .u=0 dans I (3.17)

Et l'equation(3.16) devient :

?u

+ã.ó.r = --G.ó.r (3.18)

?n

Les conditions aux limites de notre problème sont donc :

surF1 :u=-- 2.ð^(3.19)

I0

surF0UF2UF:u=0 (3.20)

?u

sur Fi U Fr: + ã.ó.r = --G.ó.r (3.21)

?n

Ainsi le problème est bien posé et on pourra appliquer une méthode de résolution numérique. Pour cette résolution, on sera obligé de remplacer les frontières à l'infini par des frontières ayant des coordonnées finies assez grandes par rapport à la longueur du casing L. La méthode de résolution utilisée est la méthode des éléments finis.

Avant d'appliquer cette méthode de résolution, une approche par perturbations a été appliquée au problème qui a permis de considérer le casing comme une frontière (ceci en réduisant son épaisseur) et donc de créer de nouvelles conditions aux limites sur le casing, cette nouvelle approche a été appliquée durant le stage de M. El Moussawi, pour plus de détails sur cette approche veuillez consulter la référence [1].

Avec les notations de la figure (3.1)

_r

S

1

a 1 au a 1 au ù.u

S_i

S₀

S

2

FIG. 3.2 notations dans le domaine I En appliquant cette méthode le problème devient

sur I, a ^{aa a}

( .)+ ( . ) -- i. .u = 0

ró1.rr z ó1.rz r

au

sur Fr, = 0,

az

au

sur S1, -- A.u = --B

ar

au

sur Si, --C.u=0

ar

au

sur S0, --A.u=0

ar

sur S2?F, u=0

A,B et C sont les constantes suivantes :

Où r¹ : le rayon intérieur du casing

r2 : le rayon éxterieur du casing

a0 : la résistivité du casing

a1 : la résistivité de la formation

U0 : l'intensité du courant injecté à travers le casing

3.2 Résolution par la méthode des éléments fi-

nis

3.2.1 Formulation variationnelle

On a donc le problème (PV) à résoudre par la méthode des éléments finis.

Soit U = {v E H¹(Ù); v = 0 sur S2 U F} l'espace des fonctions admis-

sibles pour ce problème, (PV) est équivalent à la formulation variationelle suivante

(

où (,) représente le produit scalaire dans L²(Ù).

On applique les formules de GREEN on obtient :

Or⁽¹

O a1.r^.Ou Or), v) +(O Oz⁽¹

a1.r.OuOz ),v)-(i. W.u _r.u,v)=0 (3.22)

⁽ⁱ )

a1.r^.Ou Oz, ^Ov ) + ₍1 du

Oz ) ⁰

.u, _v .
ra1.rdn

(1 .

a1.r

OuOr ,

(1 .

a1.r

dn, v)r =0 (3.23)

du

En substituant on obtient :

Alors,

( 1

ó1.r2(A.1S1?S0(z) + C.1S_i(z))u, v)r=r² = ₍B

ó1.r2 .1S₁(z), v)r=r² (3.24)

Donc le problème précédent admet pour formulation variationnelle

½ trouver u E U véri~ant

( PV) a(u,v-u)=l(v-u) pour toutvEU

Où a et l sont respectivement la forme bilinéaire symétrique et la forme linéaire données par :

?v ?r ) +(¹

ó1.r .^?u_?z,^?v

?z) + (i. ù.u _r.u,v)

a(u,v) = ₍1

?u ?r ,

.

ó1.r

+( 1

ó1.r2 (A.¹S1?S0(z) +C.1Si(z))u, ^v)r=r²

l(v) = ( B 2 .1S_i(z), v)S1

.r

^ó1

3.2.2 Forme discrete

La méthode consiste à chercher une approximation uh de u dans un sous espace Uh de dimension finie mr * mz dont on désignera par ?i,j(i = 1, ..., mr et j = 1, ..., mz) des fonctions de bases.

On peut donc décomposer uh et vh dans la base des _?i,j(r, ^z) :

u(r, z) = _Xnr _Xnz

i=1 j=1

v(r,z) = _Xnr _Xnz

p=1 q=1

ui,j?i,j(r, ^z)

vp,q?p,q(r, ^z)

Le problème discret associé à (PV) est donc
½ trouver uh E Uh vérifiant

(PV)h a(uh, vh - uh) = l(vh - uh) pour toutv E Uh

En injectant dans (PVh) les expressions de uh et de vh, il vient, en

notant wi,j⁼ vi,j^- ui,j

Les fonctions de bases étant connues, les valeurs prises par a et l sur celles-ci sont des coefficients connus. Notons les par Ai,j,p,q et par Fp,q respectivement :

Ai,j,p,q = a(?i,j, ?p,q) Fp,q = l(?p,q) i,p E (1, ..., r) j, q E (1, ..., z)

Le problème approché consiste donc à trouver r * z nombres u solutions des équations linéaires suivantes :

_Xnr _Xnz wi,j _Xnr _Xnz Ai,j,p,qup,q = _Xnr _Xnz wi,jFi,j V ù E <^nr * <nz i=1 j=1 p=1 q=1 i=1 j=1

En optant pour une numérotation des noeuds par ligne :

u 1 = u1,1, u 2 = u2,1, ..., u nr = u1,2, u nr+1 = u2,2, ..., u nr*nz = unr,nz

Les équations linéaire précédente s'écrivent sous forme condensée

trouver le vecteur U E ^<nr*nz tel que AU = F

Où la matrice de rigidité (stiffness matrix) A est de taille ( r * z) * ( r * z) (symétrique définie positive) de coefficients Ai,j,p,q , et le second membre F est un vecteur de composante fp,q
·

3.2.3 Discrétisation Q1

On réalise une rectangulation du domaine I comportant r * z noeuds et _e éléments rectangulaires :

Ù = U^neÙe

Notons k le pas de discrétisation en z et h le pas de discrétisation en r (supposés ici constants).

Les fonctions de bases ?i,j sont choisies affines par morceaux pour chaque coordonnée séparément (élément Q1) valant 1 au noeud numéro (i,j) et 0 aux autres :

?i,j(r, ^z) = ?r i (r)?^z _j(z)

(r - ri_1)

_hsi r E [ri_1,ri]

(ri+1 - r)

_hsi rE [ri, ri+1]

(z - zj_1)

_ksi z E _{[zj_1,} zj]

(zj+1 - z)

_ksi z E [zj,zj+1]

Avec h = ri - ri_1 V i = 1, .., nr et k = zj - zj_1 V j = 1, .., nz

Donc la matrice de rigidité A est construite de la facon suivante :

Ai,j,p,q = ( 1

w.u

.Vr?i,j, Vr?p,q) + (1 .Vz?i,j,Vz?p,q) + (i. .?i,j, ?p,q) -

ó1.ró1.rr

( 1

ó1 .r2 (A.¹S1?S0 (zj) + C.¹S% (zj))?1,j, ?1,q)r=r²

= ( 1 .Vr?r _i , Vr?^r _p)(?^z j, ?^z _q) + (Vz?^z _j,Vz?^z q)( 1

ó1.r.?r ⁱ , ?^r _p)

ó1.r

Pour w = 0 on a

Ai,j,p,q = ( 1 .Vr(pr _i , Vr(p^r_p)((p^z _j,(p^z _q) + (Vz(p^z _j,Vz(p^z q)( 1

ó1.r.(pri, (p^r_p)

ó1.r

-(¹

ó1.r2(A.1S1uS0(zj) + C.1S_%(zj))(p^r ₁(p^z _j,(p^r 1(pz q)r=r²

Le ehoix des fonetions de bases fait que A est pentadiagonale. Considérons les matrices T, F, M^j de taille nr * nr définies par :

Mj = [m^jt,e], aveemjt,e= ⁽¹

ó1.r2 (A.¹S1uS0(zj) + C.¹S%(zj))ä1,e, seul

m^j1,1 est non nulle.

T = [tj,q] = ((pj, (pq), done tj,q = käj,q

D = [di,p] = ( 1 1

.(pi, (pp), donedi,p = äi,p

ó1.ró1.ri

F = [fi,p] = (¹.Vr(pi,Vr(pp), done ó1.r

+

1

ó1

₍1

rt+1

1₍ 1 1

+ ) sii=p+1

ó1 rt rt+1

1₍1 1

+ )

ó1 rt rt_1

sii = p-1

1

+ ) sii=p

rt rt_1

2

-

-

1

fi,p = 2

?

? ???????????

????????????

Et B la matriee de taille nz * nz définies par

B = [bj,q] = (Vz(pj, Vz(pq), done

k si j = q avec i=61et p=61

2

k sij=q=1

1

k sij = q+1

1

k sij = q-1

1

Avec ces notations, la matrice A s'écrit par exemple avec mr = 5 et mz = 3:

(F+M1). k 2 + 1 _k.D

- 1 _k.D

0

- 1 _k.D
k.(F+M2)+ ² _k.D

-

?

? ? ? ? ? ? ? ? ? ? ? ? ?

1 _k.D0

- 1 _k.D
k.(F+M3)+ ² _k.D

Les 0 représentent des matrices nulles de taille mr * mr.

Pour le second membre, si on note V ^j le vecteur à mr éléments seul le premier coefficient est non nul et vaut ó1.r2 .¹S1 (zj), Alors F s'écrit :

B

?

? ?

? V 1_.^k

? 2

F = ? V².k

V ³.k

Nous utiliserons un gradient conjugué pour résoudre le systéme linéaire A.U = F.

La discrétisation et le gradient conjugué sont codés en Scilab sous le nom Cwmlat - bis, ce nom est inspiré du nom Cwmlat qui est le nom d'un code éléments finis chez Schlumberger.

Cwnlat résout le même problème en se basant sur la même formulation variationnelle et la même discrétisation, seulement le maillage utilisé est à pas variable contrairement à celui utilisé dans Cwnlat - bis, pour plus de détails sur le code Cwnlat consulter la référence [31.

Données numériques utilisées dans Cwnlat - bis

r1 = 0.078486 m Rayon intérieur du casing.

r2 = 0.0889 m Rayonextérieurducasing.

L nes de courant c culées pa

rinf 2540 m L'infini en r.

0

ele1 -3000 m Extrémité supérieure de l'électrode d'injection.

00

ele0 -3500 m Extrémité inférieure de l'électrode d'injection.

shoe -4500 m Extrémite inférieure du casing.

0

zinf -5300 m L'infini en z.

ó1 = 0.01 Ù Résistivité de la formation.

E 0.0001 Coefficient de régularisation.

00

FIG. 3.3 - Lignes de courant calculées par Cwnlat-bis dans une formation homogène à 0.01Ù

Chapitre 4

Méthode de Factorisation

4.1 Motivation

Au cours d'un log simulé le CHFR effectue plusieurs stations en remontant. A chaque station un courant est injecté à travers le casing, donc pour chaque station le programme Cwnlat doit résoudre sur tout le domaine Ù (qui est assez grand) le même problème (PV). En général le nombre de stations est assez élevé, c'est pourquoi le temps de calcul du programme Cwnlat pour un log simulé est très long.

A maillage fixé on se propose d'effectuer les calculs sur un sous domaine fixe Ù1 qui contient l'ensemble des positions de la sonde. Le nombre des inconnues est ainsi réduit. Le problème est alors de déterminer des conditions aux limites sur le bord t(^Ù1 est la partie du domaine qu'on souhaite enlever) qui résument exactement le comportement de la solution sur Ù1, ou autrement dit d'éliminer les inconnues sur Ù1.

Les conditions aux limites transparentes seront donc calculées sur le domaine _t, en utilisant l'opérateur Dirichlet-Neumann :

=Q.u+q. (4.1)

Où Q est un opérateur autoadjoint et q un vecteur. Ainsi grace à cette opérateur, si on suppose connaitre une solution u de la famille des solutions de notre problème sur _t de passer à une condition aux limites de type Neumann sur la frontière transparente,pour plus de détails sur cette méthode veuillez consulter la référence [21.

Ces conditions aux limites transparentes une fois calculées seront inté-

grées dans le programme Cwnlat - bis qui se base sur le cadre mathématique décrit dans le chapitre 3, et plusieurs tests seront effectués par la suite.

Dans un premier temps je présente le calcul formel de la méthode de factorisation, puis le calcul discrét avant de conclure par différents tests. La figure ci dessous donne les notations pour le développement de la méthode.

_{r s}

Ù₁

_z

S

1

z_t

t

S_i

S₀

S

2

Ù = Ù Ù

-

2

1

FIG. 4.1 Notations pour le développement de la méthode

4.2 Présentation du calcul formel de la méthode

de factorisation

? 1 ?u ? 1 ?u ù.u

On rappelle le problème à factoriser :

sur Ù, ? ?? ?

( . ) + ( . ) -- i. .u = 0

ró1.rr z ó1.rz r

?u

sur Fr(électrode de recueil), = 0

?z

?u 1

surFs, + .u=0

?z E

?u

sur S1, -- A.u = --B

?r

?u

sur Si, --C.u=0

?r

?u

sur S0, --A.u=0

?r

sur S2?F, u=0

Où A,B et C sont des constantes :

Avec r¹ : le rayon intérieur du casing

r2 : le rayon extérieur du casing

ó0 : la résistivité du casing

ó1 : la résistivité de la formation

?u 1

+

?z

U0 : l'intensité du courant injecté à travers le casing

.u = 0, le coefficient E est choisi très petit (10^-6).

E

Avant de commencer les calculs, j'aimerais insister sur le fait que l'opérateur Q et la fonction q dépendent de z, puisqu'on a choisi de réduire le domaine parallèlement à la surface. Et nous allons faire jouer à z un rOle particulier comparable à celui du temps pour le contrOle de problèmes d'évolut ion.

Soit U = {v E H¹(Fzt);v = 0 sur F} l'espace des fonctions admis-

sibles pour ce problème.

Remarque Pour les profondeurs z situées en dessous du pied de casing

on prendrait U = {v E H¹(Fz_t); v = 0 sur S2 U F}

Soit donc ? E U,on a

On dérive cette équation par rapport à la variable z :

( O Oz ( 1

ó1.r^.Ou Oz ), ?) = (OQ Oz.u, ?) +(Q.Ou Oz , ?) +(^Oq

Oz ,?) (4.3)

De (PV) on a

(i^ù.u

r

u, ?) - ( O Or ( 1Oz .u, ?) + (Q.ó1r.Q.u, ?) +

ó1.r^.Ou Or ), ?) = ₍OQ

(Q.ó1r.q, ?) +(^Oq

Oz ,?) (4.5)

On intègre par partie le terme (O Or^{( 1}

ó1.r.OuOr) par rapport à la va-

riable r, ce qui nous donne :

ó1.r2 1S1(z).?(r2) = ₍OQ

B Oz .u, ?) + (Q.ó1r.Q.u, ?) + (Q.ó1r.q, ?) + ₍Oq

Oz , ?)

Cette égalité fait intervenir les différentes conditions aux limites de notre problème. Le calcul précédent est valide pour tous les u de la famille vérifiant l'équation aux dérivées partielles de notre problème et la condition limite sur Fr. Ainsi u(r, z) est quelconque dans la famille des solutions, on

en déduit de l'identité précédente :

Et, (q,?) + (Q.ó1r.q,?) + ó1.r21S1(z).?(r2) = 0 (4.7)

B

Avec ? et ø dans U

Il nous reste à déterminer Q et q pour z = 0 (en surface). sur F₅

?u

on a la condition aux limites +

?z

1 ?u

.u = 0, et sur Fr on a = 0,

E ?z

alors d'après la relation Dirichlet-Neumann on a sur F₅ (Q + D.¹)u+ q = 0,

E

comme u est arbitraire alors on a

{ -¹_E.D surF₅

Q(0)= ₀ surFr
q(0) = 0 sur F₅

1

äi,j,

Ici D est une matrice diagonale de taille _nr'3*nr'3 telle que [di,j] ⁼ó1.ri

nr'₃ est le nombre de noeuds sur F₅.

L'opérateur Q qui dépend de z et agit sur des fonctions de la variable r vérifie une équation de Riccati qui a la particularité de ne pas avoir de terme linéaire. En réunissant les équations (4.6), (4.7) et les conditions initiales pour Q et q, on obtient le système

?

? ??????????????

???????????????

ù.u (i. _r

ø, ?)+( 1 .Vrø,Vr?) - (Q.ó1r.Qø, ?) - (?Q ?z .ø, ?)+ ó1.r

ó1.r2(A.1S1?S0(z) + C.1S_i(z))ø(r²)?(r²) = 0

1

(q,?) + (Q.ó1r.q,?) + B

1 .

ó1 .r

?u + Qu = -q ?z

ó1.r21S1(z).?(r2) = 0

1

Q(0)= -.D

E

q(0) = 0

Ce système est découplé car on peut intégrer les deux premières equations en z de 0 à zinf pour obtenir Q et q, puis u est obtenu par intégration retrograde de la troisième equation.

4.3 Forme discrete

= Q.u+qpar

1

On approche la relation Dirichlet-Neumann:

?u .

ó1.r?z

un schema differences finies décentré à gauche :

1

Ici D est une matrice diagonale de taille nr * nr [di,j] = .äi,j.

ó1.ri

Et hj+1 = zj+1 - zj.

Ona

uj+1 - _uj

(D. , ?) = (Qj+1.u^j+1 , ?) + (^qj+1, ?) (4.9)

hj+1

uj - _uj_1

(D. ,?)=(Qj.uj,?)+(qj,?) (4.10)

hj

Après une soustraction entre l'equation(4.9) et l'equation (4.10) on obtient :

hj.u j+1 - u^j(hj+1 + hj) + hj+1.u^j_1

(D., ?) = (Qj+1.u^j+1 - Qj.u^j, ?) + hj+1.hj

(^qj+1 -qj,?) (4.11)

Pour que le terme de gauche de l'equation (4.10) soit bien symé-

2

trique nous allons multiplier cette equation par :

hj+1+hj

:

(D. hj+1+hj uj+1 - 2u^j + 2hj+1

2.hj hj+1+hj _uj_1 2

hj+1hjhj+1 + hj

, ?)= [(Qj+1u^j+1 - Qju^j,?) +

(^qj+1 - qj, ?)] (4.12)

hj

hj+1+hj 2.hj

Pour alléger l'écriture posons hj+1 ₂ = 2, ^áj=hj⁺¹ +

áj.u j+1 - 2.u^j + áj.u^j-1 1

(D., ?) = [(Qj+1.u j+1 - Qj.u^j, ?)

hj+1.hj h_{j +}1

2

+(^qj+1 - qj, ?)] (4.13)

Une semi-discrétisation de l'équation aux dérivées partielle (3.17) suivant la variable z nous donne :

Où Vr désigne l'opérateur _?r, alors l'équation (4.13) devient

?

Qu'on peut encore écrire

+(Qj+1(u^j+1 - u^j), ?) + (qj+1 - qj, ?)] (4.16)

De l'équation (4.8) on a

uj+1 - u^j = hj+1(D - hj+1Qj+1)-1(Qj+1u ^j +^qj+1) (4.17)

Alors l'équation (4.16) devient

1

Une intégration par parties sur le terme ((Vr. .Vr)u^j, _?)zt par

ó1.r

rapport à la variable r nous donne l'équation suivante :

1

= [((Qj+1 - Qj)u^j, ?) + hj+1(Qj+1(D - hj+1Qj+1)^-1Qj+1u^j, ?)
hj+ ¹

2

+ hj+1(Qj+1(D - hj+1Qj+1)^-1qj+1, ?) + (^qj+1 - qj, ?)] (4.19)

On obtient donc l'équation :

+ ó1.r2 (A.¹S1?S0(zj) + C.1S_%(zj))u^j(r²)?(r²)

1

1

= [((Qj+1 - Qj)u^j, ?) + hj+1(Qj+1(D - hj+1Qj+1)^-1Qj+1u^j, ?)
hj+ ¹

2

+ hj+1(Qj+1(D - hj+1Qj+1)^-1qj+1, ?) + (^qj+1 - qj, ?)] (4.20)

L'équation (4.20) qui fait intervenir les différentes conditions aux limites de notre problème. Le calcul précédent est valide pour tous les u^j de la famille vérifiant l'équation aux dérivées partielles de notre problème et la condition limite sur z = 0. Ainsi u^j(r, zj) est quelconque on en déduit de l'identité précédente :

ù.u

hj+ 1 ₂(i.ø, ?) + h_j+ 12( 1

ó1.r.Vrø,Vr?) - (Qj+1ø, ?) + (Qjø, ?)

r

1

+hj+ 1 2 . ó.r2(A.1S1?S0(zj) + C.1S%(zj))ø(r2)?(r2)

- hj+1(Qj+1(D - hj+1Qj+1)^-1Qj+1ø, ?) = 0 (4.21)

Avec ? et ø dans U Et

B

(^qj+1 -qj,?)+hj+1(Qj+1(D-hj+1Qj+1)-1qj+1,?)+hj+1 ₂ . ó1.r2 1S₁(zj).?(r²) = 0

(4.22)

Cherchons maintenant une approximation øh de ø dans un sous espace Uh de dimension finie nr dont on désignera par îi (i = 1, ..., nr) les fonctions

de bases, choisies affines par morceaux valant 1 au noeud numéro i en r et 0 aux autres.

Tout élément Ph de Uh est donc de la forme

Où les Pri sont des coefficients (projections de vh sur les fonctions de bases), On a donc

ù.u

hj+ 1 ₂ (i.îl, îe) + hj+12 ( 1

_ó1.r.Vrîl,Vrîe) - (Qj+1îl, îe) + (Qjîl, îe)

r

1

-hj+1(Qj+1(D - hj+1Qj+1)^-1Qj+1îl, îe) +h_j+ 1 2 ._ó.r2(^A.1S1?S0(^zj)

+C.¹Si(zj))î1(r²)81,lî1(r²)81,e = 0 (4.23)

Et

B

-hj+ 1 2 .ó.r2 ¹S1 (zj).î1(r²)81,e = 0

Rappellons que Qj est un opérateur autoadjoint et que les fonctions de bases îl sont des fonctions de r, alors _Qjîl est aussi une fonction de r en effet puisqu'on applique à la fonction de bases îl l'opérateur Qj.

D'aprés le théorème Noyau de Schwartz il existe un noyau gj

défini sur zt X zt tel que

_Z

Qjîl(r) = gj(r,x)îl(x)dx

z_t

alors

_ZZ

(Qjîl, îe) = gj(r,x)îl(x)îe(r)dx.dr

supp(îl)supp(îe)

Pour réaliser la condensation de masse nous allons utiliser la méthode de quadrature pour le calcul les intégrales. Alors,

kt+1 + kt ke+1 + ke

(Qjît, îe)'₂ .₂^gj(rt,xe).

Posonsgj t,e = ke+1 + ke

₂^gj(rt,xe) (Qjît, îe) = kt+1 + kt

2gj t,e

Nous allons donc considérer comme inconnue de l'équation (4.23) la

matrice Qj = [gj t,e]^.

On obtient donc

T. Q j+1 - T . Q _j + h j+1T . Q j+1( D - hj+1 Qj+1) ^_1 Qj+1

-hj+ 1 2 (E+F+M^j) = 0 (4.24)

Et

T.^qj+1 - T.qj + hj+1 Qj+1(D - hj+1 Qj+1)_¹.T.^qj+1 - h_j+ 1 ₂ .V ^j = 0 (4.25)

Où V ^j est un vecteur de nr éléments, seul le premier élément est non nul et vaut _{ó1.r21S1(zj),} tous les autres éléments sont nuls, quant à T, F,

B

M^j, E, ce sont des matrices de taille nr * nr :

k1

T = [tt,e], avec tt,e= kt+1 + kt

2 ät,e,si ^l = e = 1 t1,1 = 2

E = [et,o], avec et,o= ^ù.u ät,o

^r

t

F = [ft,e], avec

-

-

1₍ 1 1

+ )sil = e+1

ó1 rt rt+1

1₍ 1 1

+ )sil = e-1

ó1 rt rt_1

1

+ )sil = e

rt rt_1

2

1

ft,e = 2

?

? ???????????

????????????

En ce qui concerne la matrice M^j seul le coefficient _mj1,1 est non

nul et vaut ó1.r2 (A.¹S1uS0(zj) + C.1S_i(zj)), tous les autres coefficients sont 1

nuls.

Le calcul de l'opérateur Qj et du résidu qj se fait suivant les indices croissant de zj, il faut donc exprimer Qj+1 en fonction de Qj et qj+1 en fonction de qj.

Donc en revenant aux équations (4.24) et (4.25) on a

Qj+1(D -- hj+1 Q j+1) ^-1((D -- h j+1 Q j+1) + h j+1 Qj+1)
= Qj +hj+ 2 .T ^-1(E+F+M^j)
? Qj+1(D^--hj+1 Qj+1)^-1D = Qj+ hj+ 2 .T ^-1(E + F + M^j)

? ( Q^-1

j+1 -- hj+1D^-1)^-1 = Qj+ hj+ 2 .T ^-1(E + F + M^j)

? Qj+1 = (hj+1D^-1 + ( Qj+ hj+ 2 .T ^-1(E + F + M^j))^-1)^-1 (4.26)

Et

(I + hj+1 Qj+1(D -- hj+1 Qj+1)^-1)qj+1 = qj + hj+ 2 .T -1.V ^j

? ^qj+1 = (I + hj+1 Qj+1(D -- hj+1 Qj+1)-1)-1(qj +hj+ 2 .T -1.V ^j) (4.27)

Nous venons donc de construire un schéma numérique sur l'opérateur Qj et le résidu qj, il nous reste à déterminer Q1 et q1.

Par identification avec (4.8) écrit pour j = 0 en utilisant le caractère arbitraire de u¹ on a

q1 = 0et (Q1ø, ?) = --(¹ ø, _?)8 donc

²

01

u D0 Q1⁼₀₁t)

--D¹.¹

²

Ici D¹ est une matrice diagonale de taille (mr -- m)*(mr -- m) telle que

D0 est une matrice nulle de taille m * m, où m est le nombre de noeuds de l'emplacement de l'électrode de recueil.

Et 0¹ une matrices nulle de taille m * (mr - m).

Maintenant que Q1 et q1 sont déterminés et si le courant injecté est continu alors l'algorithme de calcul de l'opérateur Q3 et le résidu q3 est

de j = 1 à mzt faire

calcul de la matrice ^M3

Q3+1 = (h3+1D^-1 + ( Q3+ h₃₊ 1 2 T ^-1(F + M³))^-1)^-1

calcul du vecteur ^V3

^q3+1 = (I + h3+1 Q3+1(D - h3+1 Q3+1)-1)-1(q3 +h₃₊ 1 ₂ .T -1.V ³⁾

fin de la boucle sur j

L'algorithme (ALGz_t) est programmé en Scilab avec les mêmes données numériques utilisées dans le programme Cwmlat - bis. L'opérateur Q et q sont sauvegardé dans des fichiers différents et seront intégré dans

Cwmlat - bis.

Une façon de vérifier le calcul de l'opérateur Q et le résidu q est de les calculer sur tout le domaine Ù suivant les indices croissants en z3 et en calculant la solution grace à la discrétisation de la relation Dirichlet-Neumann

suivant les indices décroissants en z3 u³ = u³⁺¹ - h3+1D^-1 (Q3+1u³⁺¹ + ^q3+1),

avec unzinf = 0 (condition aux limites de type Dirichlet sur F). Pour résumer

on a l'algorithme :

(^ALGzinf)

?

? ???????????????

????????????????

de j = 1 à _nzinf faire calcul de la matrice ^M3

Q3+1 = (h3+1D^-1 + ( Q3+ h₃₊ 1 ₂ T ^-1(F + M³))^-1)^-1

calcul du vecteur ^V3

^q3+1 = (I + h3+1

Q3+1(D - h3+1 Q3+1)-1)-1(q3 + h₃₊ 1 ₂ .T -1.V ³⁾

fin de la boucle sur j

de i = nzinf à 1 faire

u -1 = u - h D-1( Q u + q )

fin de la boucle sur i

L'algorithme _(ALGzinf) est aussi programmé en Scilab sous le nom cal - Q - q - zinf avec les mêmes données numériques utilisées dans le programme Cwnlat - bis. En comparant cette solution avec celle trouvée avec le programme Cwnlat - bis on a trouvé une erreur moyenne quadratique de l'ordre de 2.23 10^-5 entre les deux solutions.

La figure ci-dessous représente les lignes de courant calculées avec Cwnlatbis et cal - Q - q - zinf. Les deux autre courbes représentes la distribution spatiale des erreurs absolues et relative entre les deux solutions.

La distribution spatiale des erreurs absolues et relatives montrent une variation de l'erreur au voisinage des électrodes d'injection et de recueil, ceci est dû au fort gradient au voisinage des électrodes, on remarque aussi que la distribution des erreurs relatives montre une variation de l'erreur relative au voisinage de la frontière S2, ceci peut étre dû à une mal prise en compte de la condition aux limites sur cette frontière, dans le schéma numérique.

^{Distribution spatiale des erreurs abs Distr bution spatiale des er eu s ela}

-¹

^-5000

FIG. 4.2 Test de comparaison entre les lignes de courant calculées par

Cwnlat - bis et les lignes de courant calculées par cal - Q - q - zinf dans

une formation résistive à 0.01Ù.

4.4 Prise en compte des conditions aux limites

transparentes dans Cwnlat-bis

En réduisant le domaine Ù le problème devient :

sur ^p2= ^p - p1, ^-or⁽¹

o ó1.r^.ou or^)-o oz⁽¹

ó1.r.ouoz)+i. ù.u _r.u=0

ou oz

= Q.u+q

1

sur Fz_t(frontière transparente), .

ó1.r

ou

sur S1, - A.u = -B
or

ou

sur Si, -C.u=0

or

ou

sur S0, - A.u = 0
or

sur S2UF, u=0

4.4.1 Formulation variationnelle

Soit réduit = {v E H¹(^p₂); v = 0 sur S2 U F} l'espace des fonc-

tions admissibles pour ce problème et u, v deux éléments de réduit , Alors

-( o or( ¹

ó1.r^.ou or ),v)-( o oz ( 1

ó1.r.ouoz ),v)+(i. ù.u _r.u,v)=0 (4.28)

où (,) représente le produit scalaire dans L²(^p₂).

On applique la formule de GREEN, on obtient :

ov or ) +(¹

ó1.r^.ouoz^, ovoz) - (i.^ù.u_r.u, v)+( ¹

ó1.r^.du

(1 .

ó1.r

ou or ,

dn , v)?Ù₂^{= 0}

(1 .

ó1.r

_dn,v)z=0⁼0 (4.29)

du

En substituant on obtient :

Alors,

(1 .

ó1.r

?u ?r ,

?v ù./1

+

?r ) +(¹

r

ó1.r.?u?z, ^?v

?z) + ^{(i u v)} + ^{(Q u v)}

( 1

ó1.r2(A.1S1?S0(z) + C.1S%(z))u, ^v)r=r² = --(q, v)z=0

+(^B

ó1.r2 .¹S1 ^{(z), v})r=r² (4.30)

Donc le problème précédent admet pour formulation variationnelle

(PV)réduit { trouver u E ^Uréduit vérifiant

a(u, v -- u) = l(v -- u) pour tout v E _Uréduit

Où a et l sont respectivement la forme bilinéaire symétrique et la forme linéaire données par :

?v ù./1

)

?r ) +(¹

ó1.r.?u?z, ^?v

?z) + ⁽ⁱ⁾ + ^(Q

r

+( 1

ó1.r2 (A.¹S1?S0(z) + C.1S%(z))u, v)r=r²

a(u,v) = ₍1

?u ?r ,

.

ó1.r

l(v) = --(^q,v)z=0 + ₍B

_ó1.r2^.1S1(^z)^,v)^S1

4.4.2 Forme discrete

On utilise la même méthode utilisée dans le chapitre 3, c'est à dire chercher une aproximation uh de u dans un sous ensemble ^Uhréduit de dimension finie mr * mz, avec les mêmes fonctions de bases.

On peut donc décomposer uh et vh dans la base des ?i,j(r, ^z) :

u(r, z) = _Xnr _Xnz ui,j?i,j(r, z), v(r, z) = _Xnr _Xnz vp,q?p,q(r, ^z) i=1 j=1 p=1 q=1

Le problème discret associé à (PV) ^réduitest donc

(PV)hréduit { trouver uh E ^Uhréduit vérifiant

a(uh, vh - uh) = l(vh - uh) pour toutv E _Uhréduit

En injectant dans (PVh^réduit) les expressions de uh et de vh, il vient, en notant wi,j⁼ vi,j- ui,j, noeuds par ligne :

Les fonctions de bases étant connues, les valeurs prises par a et l sur celles-ci sont des coefficients connus. Notons-les par Ai,j,p,q réduitet par F _p,q ^réduitrespectivement :

Ai,j,p,q réduit = a(?i,j, ?p,q) F p,qréduit = l(?p,q) i,p E (1, ..., mr) j, q E (1, ..., mz)

Le problème approché consiste donc à trouver mr * mz nombres u solutions des équations linéaires suivante :

_Xnr _Xnz wi,j _Xnr _Xnz Ai,j,p,q ^réduitup_,q = _Xnr _Xnz wi,jFi,j^réduit v ù E <^nr * r:^nz i=1 j=1 p=1 q=1 i=1 j=1

Les équations linéaire précédente s'écrivent sous forme condensée trouver le vecteur ^Uréduit E<^nr*nztel que AréduitUréduit = _Fréduit

Où la matrise de rigidité (stiffness matrix) Aréduit , est de taille

réduit

,j,p,q

(mr *mz) * (mr * mz) (symétrique définie positive) de coefficients Ai

, et le second membre _Fréduit est un vecteur de composante fp,qréduit.

4.4.3 Discrétisation Q1

Avec le même maillage utilisé dans le chapitre 3, la matrice de rigidité ^Aréduit est construite de la facon suivante :

Ai,j,p,q^réduit = ( 1 w.u

.Vr?i,j, Vr?p,q) + ( 1 .Vz?i,j, Vz?p,q) + (i..?i,j, ?p,q) +

ó1.ró1.rr

(Q.?i,1, ?p,1)z=0 + ( 1

ó1.r2(A. ¹S1?S0 ^(zj) +C.¹Si(zj))?i,j,?p,q)r=r²

= ( 1 .Vr?r _i , Vr?^r_p)(?^z _j,?^z _q) + (Vz?^z _j, Vz?^z q)( 1

ó1.r.?r ⁱ , ?^r _p) +

ó1.r

₍1

ó1.r2 (A.¹S1?S0(zj) + C.1S_i(zj))?^r ₁?^z j, ?^r 1?z q)r=r2 +

w.u

(Q.?r i ?z _1, ?r p?^z1)z=0 + (i. .?i,j, ?p,q)

r

Pour w = 0 on a

Ai,j,p,q^réduit = ( 1 .Vr?r _i , Vr?^r _p)(?^z j, ?^z _q) + (Vz?^z _j, Vz?^z q)( 1

ó1.r.?r ⁱ , ?^r _p)

ó1.r

+(¹

ó1.r2 (A.¹S1?S0(zj) + C.1S_i(zj))?^z ₁?^z j, ?^r 1?z q)r=r2 +

(Q.?r _i , ?^z ₁, ?p?^z1)z=0

Avec les mêmes notations utilisées dans la chapitre 3 la matrice Aréduit s'écrit par exemple avec mr = 5 et mz = 3 :

1

(F+M1). k 2 + _k.D+Q*k

- 1 _k.D 0

k.(F+M2)+ ² _k.D

- 1 _k.D

- 1 _k.D

k.(F+M3)+ ² _k.D

0

- 1 _k.D

?

? ? ? ? ? ? ? ? ? ? ? ? ?

?

? ? ? ? ? ? ? ? ? ? ? ? ?

Aréduit =

Les 0 représentent des matrices nulles de taille r * r.

Pour le second membre, si on note V ^j le vecteur a _r éléments seul le

premier coefficient est non nul et vaut ó1.r2 .¹S1 (zj), Alors ^Fréduit s'écrit :

B

? ? ?

? V ¹.^k ₂+ k.q

Fréduit = ? ? V².k

V³.k

Le systéme linéaire obtenu est plus rapide à résoudre car on a diminué le nombre d'inconnues.

4.5 Tests

Pour mettre en évidence le gain en temps de calcul, on a calculé à plusieurs profondeurs des conditions aux limites transparentes qu'on a intégrées dans Cwmuat - bis et on a recupéré à chaque fois le temps que celui ci

1

met à s'exécuter. Ce temps de calcul n'est pas vraiment objectif car je n'ai pas tenu compte du caractére creux de la matrice de rigidité dans les deux

160

codes, Cwmuat - bis et Cwmuat - bis - avec - comd - uim - tram.

La courbe 4.3 représente la variation du temps de calcul de Cwmuat-bis

4

prenant en compte les conditions aux limites transparentes calculées à différentes profondeurs en fonction des profondeurs,

0

FIG. 4.3 Variation du temps de calcul de Cwmuat - bis avec des conditions aux limites transparentes calculées à plusieurs profondeurs en fonction de l'emplacement des profondeurs

Attention! Le temps t pour z = 0 représente le temps de calcul de Cwnlat - bis sans conditions aux limites transparentes.

On remarque que pour une profondeur proche de la surface (un pas en dessous de la surface) Cwnlat - bis avec conditions aux limites transparentes calculées à cette profondeur, met plus de temps à résoudre le problème que Cwnlat - bis sans conditions aux limites tranparentes. Ceci vient du fait que le nombre d'inconnues diminues peu, donc la dimension de la matrice de rigidité ^Aréduit reste importante de plus ^Aréduit est localement pleine car Q est une matrice pleine . Par contre plus on s'éloigne de la surface plus le nombre d'inconnue diminues donc la dimension de la matrice de rigidité Aréduit devient plus petite ce qui rend l'éxécution du code Cwnlat - bis avec conditions aux limites tranparentes plus rapide d'ou la décroissance en temps en fonction de la profondeur sur la courbe.

On s'est intéressé aussi aux erreurs entre les deux codes Cwnlat - bis avec conditions aux limites transparentes et Cwnlat - bis sans conditions aux limites transparentes.

Soit U(r, z) la solution calculée par Cwnlat - bis sur tout le domaine I

et U(r, z)^zt la solution calculée par Cwnlat - bis - avec - cond - aux - limi - trans sur le domaine réduit I - I1. Pour zt = {-100m, -500m, -1000m -

1500m, -2000m, -2500m}, nous avons calculé l'erreur moyenne quadratique:

U(r,z) - U(r,z)zt 2à chaque emplacement de la frontière transparente zt.

U(r,z) ₂

La figure 4.4 représente la variation de l'erreur moyenne quadratique en fonc-

tion de zt

FIG. 4.4 - Variation de l'erreur moyenne quadratique en fonction de zt. Pour savoir comment et où varie l'erreur, nous affichons pour chaque

zt les solutions U(r, z) et U(r, ^z)réduit pour zinf < z < zt, la distribution

spatiale des erreurs absolues ainsi que la distribution spatiale des erreurs relatives, sur les figures 4.5, 4.6, 4.7, 4.8, 4.9, 4.10, 4.11.

Ces tests montrent que l'erreur est systématique et varie au voisinage de l'emplacement de la frontière transparente. Ceci est dû au fait que dans le calcul de la condition aux limites transparente, la condition aux limites sur Fr est approchée par un schéma décentré à gauche, alors que dans le programme Cwnlat - bis, cette condition aux limites est approchée par un schéma centré.

^{Distribution spatiale des erreurs abs es Distribution spatiale des erreu s relat ves}

FIG. 4.5 - La condition aux limites transparente est calculée à une profondeur de -100m. Les lignes de courant calculées par Cwnlat--bis sont affichées pour

z <= --100m

^{Distribution spatiale des erreurs abs es D stribution spatiale des erreu s rela ves}

12

FIG. 4.6 - La condition aux limites transparente est calculée à une profondeur de -50Dm. Les lignes de courant calculées par Cwnlat--bis sont affichées pour

z <= --500m

Distribution spatiale des erreurs abs es x 10^-3

^{Distribution spatiale des erreu s relat ves}

FIG. 4.7 - La condition aux limites transparente est calculée à une profondeur de -1000m. Les lignes de courant calculées par Cwnlat -- bis sont affichées

pour z <= --1000m

Méthode de Factorisation

4000

^{Distribution spatiale des erreurs abs es}

x ¹⁰

Distribution spatiale des erreurs abs es x 10^-3

^{Distribution spatiale des erreu s relat ves}