dgues Et Barrages en Terre de Faible Hauteur

N.B. La durée des travaux d'injection a été de trois semaines (y compris amenée et repli) pour un coût en 1982 de 335 kF TTC) (longueur de forage «225 m, volume de coulis « 3 3 m ). Le coût de la construction de la digue a été de 1800 kF TTC (1980). 3

Fig. 9 — Principe d'imperméabilisation par injection.

18

i

Charge hydraulique Im)

Tubes crépines : K, L Piézomètres P A C : J, M , N

hydraulique dans certaines parties du profil soit supérieure à la charge dans la retenue ;

Charge dans la retenue

— soit, l'écoulement à travers le profil est alimenté, non seulement par la retenue, mais aussi par une source indépendante de la retenue, située sous la base de la digue.

1980

1981

1982

Fig. 10 — Variations de la charge hydraulique dans le profil 3 entre le début de la première mise en eau (juillet 1980) et la fin de 1982.

Les pressions interstitielles mesurées en fin de construction dans chaque piézomètre du profil 2 (elles sont comprises entre 16 et 50 kPa) étant inférieures à celles mesurées en fin de première mise en eau (c'est-àdire plus de dix mois après la fin de la construction) ou du même ordre, la première hypothèse ne peut être retenue. E n revanche, la seconde est cohérente avec les observations faites au début des travaux. Il est, en effet, extrêmement probable que la source, découverte à l'emplacement de la conduite de vidange et qui est vraisemblablement alimentée par des eaux provenant du massif constituant la rive droite, perturbe les conditions hydrauliques dans les profils voisins, notamment dans le profil 2 situé à proximité immédiate de cette conduite (à environ 4 m).

Influence des injections sur la pression interstitielle

b. l'écoulement théorique a été étudié dans l'hypothèse d'un rapport d'anisotropie de 9. Fig. 11 — Comparaison des écoulements réel et théoriques dans le profil 3 à la date du 22 novembre 1982. Les réseaux théoriques ont été déterminés en supposant que la fondation de la digue est imperméable.

Piézomètres P A C : A , B, C, D, E , F A

°

o-°

A o

Bien qu'elles aient été peu nombreuses et qu'elles aient intéressé essentiellement le profil 3, les mesures effectuées après la fin des travaux d'injection mettent en évidence l'effet bénéfique de ces derniers sur le comportement hydraulique de l'ouvrage à l'extérieur de la zone injectée (comprise entre les profils 4 et 7, fig. 5). Les mesures réalisées dans le profil 3 montrent, en effet, que les pressions interstitielles et les niveaux d'eau dans les tubes se sont abaissés pendant les deux mois suivant les travaux d'injection, alors que le niveau d'eau s'élevait dans la retenue (fig. 10). Cela implique, soit que le coulis d'injection ait colmaté des fissures s'étendant jusque sous le profil 3, soit que, avant traitement, l'écoulement dans le profil ait été alimenté partiellement par la partie de digue située à sa gauche et que les injections aient eu pour effet d'annuler ou de réduire cette alimentation latérale.

Anisotropie de perméabilité du corps de digue L a méthode la plus fiable de détermination du rapport k d'anisotropie de perméabilité — (k et k étant respec-

Charge hydraulique (m)

h

v

K

tivement les coefficients de perméabilité horizontale et verticale) du sol constitutif d'un ouvrage hydraulique en service consiste : — à mesurer les valeurs de la pression interstitielle en divers points de l'ouvrage lorsque celui-ci est le siège d'un écoulement permanent; — à calculer en ces points les valeurs de la pression interstitielle en régime permanent pour diverses valeurs de k ; h

1

1—i—i—i

—i—i—i J

A

S

O N 1980

D

J

1

1 1—l—r-T—1 1 J

A

S

O N 1981

D

1—I—V-i—r 1 — > J

A

S 0 1982

N

Fig. 12 — Variations de la charge hydraulique dans le profil 2 entre le début de la première mise en eau (juillet 1980) et la fin de 1982.

k — à retenir la valeur de —

K

assurant

la

meilleure

concordance entre les valeurs théoriques et expérimentales de la pression interstitielle. 19

L'application de cette méthode à la digue de L a Piagne se heurte à plusieurs difficultés. D'une part, aucune série de mesures de pression interstitielle ne correspond à un régime d'écoulement permanent, d'autre part, les conditions aux limites de l'écoulement sont très mal définies dans le profil 2 (existence d'une source dans la fondation à proximité du profil) et ne sont pas connues avec certitude dans le profil 3 (possibilité d'alimentation latérale ou par des fissures de la fondation). Dans ces conditions, on ne pouvait procéder qu'à une estimation assez peu précise du rapport d'anisotropie. Cette estimation a été faite à partir des mesures effectuées dans le profil 3 à la fin du mois de novembre 1981, c'est-à-dire plus de deux mois après l'achèvement des travaux d'injection, alors que le niveau de la retenue n'avait pas subi de fluctuation sensible depuis un mois. A cette époque, l'injection de la fondation avait amélioré l'étanchéité de cette dernière au droit du profil 3, et/ou réduit l'alimentation latérale de l'écoulement dans ce profil. Les résultats des mesures sont comparés sur la figure 11 aux résultats de l'étude théorique de l'écoulement permanent effectuée en supposant que le sol de fondation est imperméable et en affectant successivement les valeurs 1 et 9 au rapport d'anisotropie. O n constate que si les résultats expérimentaux ne concordent pas parfaitement avec les résultats théoriques obtenus dans l'hypothèse d'isotropie, les écarts restent relativement limités. E n revanche, les résultats expérimentaux et théoriques divergent fortement lorsque ces derniers correspondent à un rapport d'anisotropie de 9. Cela indique que le comportement hydraulique de la digue s'apparente à celui d'un ouvrage isotrope ou faiblement anisotrope *.

20

Ce résultat, qu'il convient de retenir avec une grande prudence compte tenu du peu de mesures significatives effectuées, est sensiblement différent de ceux présentés par Shérard et al. (1963)** qui recommandent d'adopter un rapport d'anisotropie de 25 pour l'étude des projets de barrages en terre construits au moyen de matériaux graveleux extraits d'emprunts homogènes.

CONCLUSION L'exemple de la digue de L a Plagne fait clairement ressortir la nécessité d'injecter le massif rocheux constituant la fondation d'une digue fermant un bassin de stockage, même si l'ouvrage est de faible importance et si une reconnaissance succincte semble indiquer que le massif est peu fissuré. L'interprétation des mesures piézométriques effectuées dans un profil de la digue semble montrer que celle-ci est isotrope ou faiblement anisotrope vis-à-vis de la perméabilité, bien qu'elle soit construite à partir d'un sol graveleux; ce résultat est toutefois à considérer avec prudence compte tenu du faible nombre de mesures significatives effectuées.

* Des réseaux d'écoulement théoriques (non représentés ici) établis pour des rapports d'anisotropie supérieurs à 9 montrent que les pressions interstitielles correspondantes sont toutes inférieures à celles obtenues pour — = 9 et, par suite, diffèrent davantage des K valeurs mesurées. * * S H E R A R D J . - L . , W O O D W A R D R . J . , GIZIENSKI S. F . , C L E V E N G E R

W. A., Earth and earth-rock dams, John Wiley and sons, Ed. (1963).

.IL

Digues et barrages en terre T^A^f/A de faible hauteur

Application de la méthode des éléments finis aux écoulements en milieu poreux Pierre HUMBERT Chef de la section modèles n u m é r i q u e s Laboratoire central des Ponts et Chaussées

RÉSUMÉ Après un rappel des équations de base régissant les écoulements en milieu poreux, une technique de résolution par la méthode des éléments finis de l'équation de Poisson non linéaire en régime transitoire est présentée. Les adaptations à apporter à cet algorithme général pour résoudre certains problèmes particuliers sont exposées et illustrées par des exemples dans le cas des écoulements en milieu poreux saturé avec surface libre, des écoulements en milieu poreux non saturé et dans le cas des calculs de nappes aquifères. MOTS CLÉS : 42 • Éléments finis (méthode) -

Écoulement (fluide) - Calcul - Nappe aquifère • Sol - Saturation - /Milieu poreux.

Si les premières applications de la méthode des éléments finis ont été orientées vers la résolution des équations d'équilibre en élasticité ou en élastoplasticité, il est très vite apparu que cette méthode pouvait aussi être utilisée pour la résolution de nombreux problèmes physiques régis par des systèmes d'équations aux dérivées partielles. E n particulier, les problèmes de répartition de la charge hydraulique dans les écoulements en milieu poreux qui, tout comme les problèmes de répartition de température ou de potentiel électrique, se ramènent à la résolution d'une équation du type A(u)=f, peuvent être traités aisément par la méthode des éléments finis. Après avoir posé les équations du problème, on abordera la technique générale de résolution et on verra sa particularisation pour trois importants types de problèmes : les écoulements en milieu poreux non saturés, les écoulements en milieu poreux saturés avec surface libre et les calculs de nappes aquifères.

21 B u l l , l i a i s o n L a b o P. et C h .

132 - j u i l . - a o ù t 1 9 8 4

Réf. 2 9 2 0

I - ÉQUATIONS DES ÉCOULEMENTS EN MILIEU POREUX

Le terme p n g ^ P y — P + - J est le coefficient d'ems

magasinement spécifique C .

1.1. Cas du milieu poreux saturé

s

L'équation générale qui régit les écoulements en milieu poreux saturé est obtenue à partir des trois équations suivantes :

L'équation générale de l'écoulement en milieu poreux saturé est donc :

di\(K

(p«) + p î = 0

(1)

f?J+ | [(1 -n) dt

par, nS p ]=0

(2)

[grad p + p g grad z]

(3)

s

— loi de Darcy : U-nU = s

avec :

les

équations (1)

n

: porosité,

q

: débit volumique,

0

: vitesse apparente du fluide,

T

U

: vitesse réelle moyenne du solide, : masse volumique du solide, : pression, : cote,

g

: accélération de la pesanteur,

k

: perméabilité intrinsèque,

p

: viscosité dynamique.

Une démarche analogue à celle du cas du milieu poreux saturé conduit à l'équation : 9 _

ÔQ\Ôh

- 7Tiir. Cs+

n

+q

— =

dt dt

dty J dt

1- Û,. grad /

Fig. 1 Courbe de perméabilité.

) ,on

)

obtient : p - (grad/j + p g grad z)

dp p = n — + — dt l-n

dn dt

P"dp — +P 1

(4)

Sans entrer dans le détail des calculs, il apparaît que, moyennant certaines hypothèses simplificatrices, l'équation (4) peut se mettre sous la forme : d i v ( X grad h) = pngl

p>-p + s

dh

et + q

avec : K= p g - : tenseur de perméabilité, M*

P

r

: compressibilité du fluide,

P

s

: compressibilité des grains solides,

a

22

(7)

L a perméabilité relative k et la teneur en eau 0 sont des fonctions de O n détermine expérimentalement la courbe K(G) (fig. 1) et la courbe v|/(0) ou courbe de succion (fig. 2). Cette dernière courbe présente

Si l'on porte dans (1) la valeur de U donnée par (3) et si l'on fait intervenir la notion de dérivée particu-

div

par

r

p z

et k

r

avec : \|/ = h — z.

p

et (3)

k .k dans (3)(k est la perméabilité relative comprise entre 0 et 1). Le coefficient S est relié à la teneur en eau 9 et à la porosité n par 0 = n S.

r

: masse volumique du fluide,

s

dans

div (k K grad h) --

p

s

I I . 2. Cas du milieu poreux non saturé Si S est le coefficient de saturation, on remplace n

— équation de continuité du solide : s

(6)

dt

et

div [p (1 -n)

Yq

s

— équation de continuité du fluide : div(p

dh

grad h) = C

: compressibilité du milieu poreux.

(5)

Fig. 2 Courbe de succion.

À

conduction électrique, électrostatique, torsion cylindrique, etc.

souvent des phénomènes d'hystérésis. Si on les néglige, on peut déterminerfe (\|/)et 0(\|i) qui sont alors introduits dans l'équation (7). r

E n régime permanent, l'équation (7) devient : div (k (v|/) K grad h) = q

(8)

r

O n retrouve le cas du milieu poreux saturé en faisant

L'équation (11), valable sur un domaine Q (fig. 3), est associée à trois types de conditions aux limites sur le contour : — condition de potentiel imposé ou condition de Dirichlet h = h(P) sur S

u

— condition de flux imposé ou condition de Neuman d'où : div (K grad h) = q qui peut s'écrire :

dxj (

" ôx

'

t

-k'Jp«, = / ; (P) sur S , dxj 2

— condition mixte, dite condition d'échange ou de Cauchy

^

-k ^-n ii

Si, de plus, il s'agit d'un milieu isotrope de perméabilité k, on obtient l'équation de Poisson : -àh=f

(9)

E n régime transitoire, si l'on pose

i

=

E(P).[h-h ]snrS , e

3

OXj

avec : h(P)

: potentiel imposé au point P,

n

: composantes du vecteur normal au contour,

{

fi (P) : flux imposé au point P, C(v|/) =

E (P) : coefficient d'échange au point P, h

e

: potentiel extérieur.

l'équation (7) s'écrit : div (k

r

K grad h) =

e«0

C + C(v|/) s

ôh ~ôt

+q

ôh Le coefficient de — comporte un terme dépendant des compressibilités du fluide, des grains solides et du milieu poreux, et un autre dépendant de la présence d'une zone non saturée, qui est souvent prépondérant par rapport au premier.

Fig.

3

E n faisant k (\|/) = 1 et pour un milieu isotrope, on obtient une équation de la forme : r

-Ah=f+c

ôh

(10)

~dt

Cette équation n'est autre que l'équation de la chaleur en régime transitoire.

1 — O n ramène le problème de la résolution d'un système différentiel à celui de la recherche du minimum d'une fonctionnelle (principe variationnel).

II - RÉSOLUTION PAR LA MÉTHODE DES ÉLÉMENTS FINIS L a méthode décrite ici est une méthode de résolution de l'équation de la chaleur : - d i v ( K grad

h)=f+c

ôh ~ô~t

qui pourra s'appliquer aux problèmes suivants : — transfert de chaleur, — écoulements en milieu poreux,

L a démarche suivie pour arriver à une formulation « éléments finis » est une démarche en deux étapes :

(H)

2 — O n ramène ensuite le problème de la détermination d'une fonction inconnue minimisant cette fonctionnelle à celui du calcul d'un nombre fini d'inconnues, par décomposition du domaine Q en un nombre fini d'éléments sur chacun desquels on impose le type de variation de la fonction inconnue. O n traitera d'abord le cas du régime permanent en milieu isotrope puis anisotrope. O n abordera ensuite le régime transitoire et la prise en compte de la nonlinéarité éventuelle. 23

II. 1. Cas du régime permanent

<">(«)= j" Z ( ^ ) d î î - 2 J ufàn

(15)

2

Principe variationnel Fonctionnelle associée à un système différentiel — condition — =f (P) sur S (Neuman); dn Su dty , on a — = — = / , ; il en resuite : dn dn t

—Ah=f

E n régime permanent, l'équation forme plus générale :

est de la

Au=f

(12)

où A est un opérateur défini, positif et symétrique agissant sur un espace de fonctions D définies sur Q. A

• Conditions aux limites homogènes

- 2

w/dn-2

Si on se limite aux fonctions de D nulles sur le contour, on démontre que, si ( 1 2 ) a une solution, elle est unique et minimise la fonctionnelle

M / Ì dS

(16)

A

(13)

F(u) = (Au,u)-2(u,f)

ôu \-E (u—u„) = 0 sur S (Cauchy);

— condition ôn

où (u, v) est un produit scalaire. O n prend en général u.vdCi.

(u, i;) =

si u rend F(u)

Réciproquement,

ôu h E u=f

elle est de la forme

(P) sur S.

2

dn du

minimum, il vérifie ( 1 2 ) .

\-Eu=

Comme

— +E du

dn • Conditions aux limites quelconques dites « non homogènes » O n montre que si \|/ est une fonction de D qui vérifie les conditions aux limites, résoudre ( 1 2 ) revient à rendre minimum

- 2

i \ Sx,

n

v|/ = /

sur S, on a :

2

f M,ifàfì Jn

.

A

(A u, « ) - 2 («, / ) + (u, A*)-

Eu

2

dS-2

uf dS

(17)

2

(*|r, A u) Généralisation à Téquation — div{K grad u)—f

E n général, cette fonctionnelle peut s'écrire :

Dans le cas d'un tenseur non isotrope a>(u) = (Au,u)-(2u,f)+

f N(u)dS

(14)

Application à Téquation — Au =f

~ [ k

K

D'après ce qui précède, la fonctionnelle O(w) est de la forme :

1

k y

2

22

on remplace l'opérateur laplacien —À par

* ( « ) = ( - A « , « ) - 2 ( u , / ) + («, -Av|/)-(v|/, - A M ) Sxj \

Sx

t

)

E n intégrant par partie (-A«,«)«-

uu, dQ, it

(u, — A\|/)

et

-Au),

O n prend alors pour fonctionnelle :

J(u) =

$ («) s'écrit :

1

Sxt Sxj

n

M u,idn-2


-

ufdCl

; i

Jn f u^dSJ <5n s

fu

Jn f „*+dS + f J 5n JJs s

5

Sn

Selon le type de conditions aux limites et sachant que i|f vérifie ces conditions, on obtient : — condition u = û~(P) sur S (Diriehlet), comme u = v|/ = u s u r S, on peut écrire

24

N (u) dS.

dCì-

Le dernier terme dépend des conditions aux limites (flux imposé ou condition d'échange) : |

N(u) dS=

j

uf dSx

|

( ^

-uf

Js

2

flux

échange

2

) dS

Si l'on pose :

Formulation matricielle Formulation locale

de; dri dÇ

Après décomposition du domaine à étudier en sousdomaines élémentaires appelés éléments (fig. 4), on approche sur ces éléments le potentiel inconnu M par un développement polynomial : u=Nl/ avec : N : matrice des l'inconnue, U

'FVE=

f

FSE=

f

£ ,

fonctions

d'interpolation de

(20)

G (N) dS

E

on peut écrire :

: vecteurs des valeurs nodales du potentiel.

E

j\àZ,àr\ dÇ

/ N |

JE ( « ) =

\ 'U

K

E

U

E

- ( ' F

E

V

+ 'F )

E

U

SE

E

Pour exprimer que J (u) est minimum, calculons SJ : E

E

6J = >U K E

E

5t/ -CF

E

£

+'F

K £

S £

) 6U

E

E n écrivant que 8J est nul pour tout bU admissible satisfaisant les conditions aux limites et en transposant l'expression obtenue, il vient : E

Fig. 4 Décomposition d u domaine.

E

K U —F E

E

VE

(21)

+ F, - SE

avec : K

: matrice de rigidité élémentaire,

E

F

: vecteur « forces de volume » élémentaire,

VE

F

: vecteur « forces de surface » élémentaire.

SE

L'équation (21) exprime la relation très simple entre le vecteur des inconnues nodales U et le vecteur des forces nodales F + F . E

VE

L a géométrie de l'élément est définie par les coordonnées des nœuds de l'élément contenues dans le vecteur X et par les fonctions de forme F qui sont les fonctions d'interpolation des coordonnées. Si r\, Ç sont des coordonnées locales, on a :

SE

Formulation globale

E

Z=[x,v,z] = F ( c ; , q , Q Z

£

O n définit les potentiels nodaux U de l'élément numéro k à partir du vecteur U des potentiels nodaux relatifs à l'ensemble des éléments par une relation de la forme : E

G

Lorsque F et N sont identiques, l'élément est dit « isoparamétrique ».

U = P" U E

(22)

G

où P représente une matrice « d'extraction » formée de 0 et de 1. De la même manière, on a : k

Si Jf et Jf sont les dérivées des fonctions d'interpolation N par rapport aux coordonnées globales et locales, on a : c

L

F =

PF

k

k

E

G

J ( u ) étant en fait l'énergie potentielle élémentaire, écrivons que la variation de l'énergie des « forces » internes de l'ensemble des éléments est égale à la variation de l'énergie des « forces » externes F pour tout 8U admissible : E

J désignant la matrice jacobienne.

G

G

S 'UÌ K

F, TI

k

L'élément de volume dCl s'écrit d£î = | J | d£, dr\ dÇ, où | J | est le déterminant de J.

E

SU* = £ «F* ÒU = 'F k

E

G

W

G

D'après (22) : X 'U X* ÒU = X 'U 'P k

k

E

E

G

k

P

k

k

bU

G

k

L a fonctionnelle s'écrit alors au niveau élémentaire :

= 'I/ (£'F*KÏ/*)5I/ 0

G

k

Jf

x

-

f fNd£l U Ja E

E

E

KJf d£i

c

! Js

c

1 U

E

R

G(N)dS U E

E

E

O n a donc : •U K 6U G

(19)

G

G

=

'F ÒU G

G

pour tout 8U admissible G

25

et par conséquent : KoU

=F

G

(23)

G

avec : k

IK^X'P*

K P k

E

k

Fig. 5

L'équation matricielle (23) représente un système linéaire dont les coefficients sont les termes de la matrice de « rigidité » globale KQ, les inconnues les composantes du vecteur des potentiels nodaux U et le second membre les composantes du vecteur des « forces » ou débits nodaux F . G

Les dérivées de N par rapport aux coordonnées locales s'écrivent :

G

r-N,u = r i o- i i

Application au triangle à interpolation linéaire T 3

|_N,r|J

O n choisit pour coordonnées locales les coordonnées barycentriques \, r\, Ç (fig. 5) : ^aire(MJK)

6

[

(

U

11 =

Ç

r

e[0, 1]

aire (IJK)

3VCC X / j f

= ^ Ï Ï K T ' aire (IJK) T|, Ç sont reliés par la relation £ + n. + Ç = 1. 6

[

0

1

1

- 1 _

L a matrice jacobienne a pour expression

]

aire (IJK) aire ( M K I )

l_0

=

Xj

T IK

yiK

"I

L

yjK

J

X

=

^K.'

JK

X

Son déterminant a pour valeur | J | = 2 S où S est l'aire de l'élément IJK. O n a d'autre part :

]

yjK

O n choisit pour fonction d'interpolation N et pour fonctions de forme F les coordonnées barycentriques elles-mêmes :

2S

1

0~

0

1

"ir*,

X

_ - y ¡ K

-1

L

JK

yjK

J\ dei dp

Jf \

1

-1

J

E

'JTL'J-'KT

2

IK

O n en déduit l'expression de la matrice élémentaire K d'après (20) :

F=N=K,n,a

4 S

X

o i r ^ -y iri

JL o

1K

X

o - i dSl

IK

k

J L-*«

J

IK

X

Lo

1

-i

J

Soit +

K

E E

= 4S

JK

x2

:

k y yKi x

'. + k x

symétrie

y

k yli x

JK

JK

x

KI

+ ky ll x

: '

k y y,j +k Xj x¡j x

J K

K

y

(24)

KyKiyu + ky X

XjJ

K I

_

symétrie

symétrie

• k yj + k x

2

y

ÌJ_

x

où e représente l'épaisseur de l'élément. Cette expression exprime directement les coefficients de la matrice élémentaire K en fonction des coordonnées des nœuds de l'élément et des coefficients du tenseur de perméabilité ou de conductibilité. Le calcul de l'intégrale sur l'élément n'est pas toujours possible E

26

analytiquement et nécessite très souvent le recours à une technique d'intégration numérique. L a figure 6 représente quelques types d'éléments utilisés dans les systèmes PPR et R O S A L I E et dans le programme CESAR.

Fig. 6. — Types d'éléments utilisés dans les PPR, ROSALIE et CÉSAR.

II. 2. Cas du régime transitoire

Après assemblage, on obtient l'équation matricielle globale :

Principe variationnel et équation matricielle E n utilisant les mêmes procédés que pour le cas du régime permanent, on montre que résoudre : u)=f+g

- d i v ( X grad

du dt

équivaut à rendre minimale la fonctionnelle : 1

y (M, t) =

2J

k« * dxi

*

KQ

+ C

G

U

G

=

(25)

F

G

Procédé d'intégration dans le temps Il existe de nombreuses méthodes d'intégration dans le temps, tant implicites qu'explicites. Parmi les méthodes implicites, on a choisi celle des différences finies dans le temps et celle de Galerkin. Méthode des différences finies dans le temps f

dxj )dS

OU A t fixé, si l'on pose Ù = — - , on obtient l'équation dt matricielle élémentaire : E

KE UE + C U E

E

{

N(M

où, comme pour le régime permanent, le dernier terme dépend des conditions aux limites (voir (18)).

=

Ul

E

f,_

At\_ 1

+

UjtQ-U&.J At

E n reportant ces valeurs dans l'équation (25), on obtient la formule de récurrence de Crank-Nicholson :

K -

F

2

où K , U et F ont la même signification que dans le cas du régime permanent, et où apparaît en plus la matrice C « d'emmagasinement », qui a pour expression : E

0

A partir des valeurs de U aux temps t _ , et t , on calcule les valeurs de U et Ù au milieu de l'intervalle de temps :

dfì

[

E

V

+

C —

At

E

+

(26)

E

C

E

=-

g ' N N dn

£

Méthode de Galerkin C'est une méthode d'éléments finis dans le temps qui nécessite le choix de fonctions d'interpolation du 27

temps dans l'intervalle [ t - i , t j . O n prendra des fonctions d'interpolation linéaires définies par : f

N

où X =

=X

tt

O n en déduit

Les valeurs interpolées des vecteurs U et F dans l'intervalle de temps sont données par :

Û(t) =

N _ U(t _ )+N U(t ) u

l

i

1

ti

l

(27)

F(t) = N , _ F ( t _ ) + i V F ( 0 |

1

(

1

( |

L'application du principe variationnel revient à écrire que sur tout intervalle t j , on a :

f

6

1 F >-i • * * 3

f

r

51/,. = 0,

quel que soit 5(7;. E n regroupant au premier membre les termes en l/j, on obtient au niveau global après assemblage la formule de récurrence de Galerkin :

L3

2 Al

J

l 5

5
(«, t) =

^-i+2F,

0

pour tout 5u(t) admissible de

(28)

à t. ;

Technique de calcul Les deux méthodes décrites ci-dessus se ramènent toutes deux à une formule de récurrence de la forme :

Or on a : 5
E

E

E

E

E

E

oU , E

étant calculé à partir de (27) :

E

BU(t)

Les matrices K et K sont fonctions de At. L'algorithme de calcul est décrit sur la figure 7. O n constate que à chaque modification de At, il faut recalculer ou pour le moins réassembler les matrices K et K et triangulariser K Les deux méthodes donnent de très t

BN . .

BN,.

t

2

t

2

v

Posons : U(td=U :

L / ( r _ ) = t/,-_ et i

1

ù = -u _ E

{

i

t

1

1

+

u ^

= Bt

l

^i

At

D'autre part, 8U est interpolé par 5t7 = ( l - > . ) 8 [ 7 _ + X,S[/

Calcul de K

E

£

(

1

État initial V (to) t

i

Calcul du second membre Parmi tous les 3U admissibles, choisissons ceux pour lesquels : E

9

§[/,._! =0

=>

dU = X6Ui E

Non

Oui

O n peut alors écrire : SO (u, t) =

Jt,' ¡ - 1

At

[(1 — À.) 'U _ i

+ At

o

• At

I

1

+ 'k 'C7J K X SfJi dX

C X SfJj àX

At 1

[(1 -X) Fi.^X

F J X oUi dX

At

1+1

Jo

X(l-X)dX=

Or, o 28

1 - et 6

X dX= o 3

Non

= Af, 9

Oui

2

Fig. 7. — Algorithme de calcul en régime transitoire (modèle linéaire).

bons résultats pour les petits pas de temps et sont inconditionnellement stables. Toutefois, pour des valeurs croissantes des pas de temps, la méthode de Galerkin oscille moins que celle de Crank-Nicholson. C'est pour cette raison que nous l'avons retenue pour les groupes 03 et 04 du système R O S A L I E et pour les programmes B I D I F et N S A T du système PPR.

Prise en compte de la non-linéarité de comportement

— la prise en compte de caractéristiques variables en chaque nœud rend obligatoire l'utilisation de la technique d'intégration numérique pour le calcul des matrices élémentaires K et C , même si les éléments sont géométriquement « droits ». E

E

L'utilisation de méthodes explicites, pour lesquelles il n'est pas nécessaire de faire une résolution à chaque pas de temps, supprimerait le premier inconvénient, mais ces méthodes ne sont pas inconditionnellement stables et nécessitent l'emploi de très nombreux pas de temps.

Présentation du problème Technique de calcul Que l'on soit en régime permanent ou transitoire, les méthodes présentées aux paragraphes II. 1 et II. 2 de ce chapitre supposent que les caractéristiques des matériaux ne varient pas, en particulier en fonction de la valeur du potentiel (charge, température) et éventuellement du temps. S'il n'en est pas ainsi, les matrices K et C dépendent de U et de t, et il est alors nécessaire de faire un calcul itératif (à chaque pas de temps en régime transitoire) en réactualisant ces matrices, jusqu'à stabilisation.

L'algorithme de calcul est décrit sur la figure 8. Par rapport au cas linéaire, on constate la présence d'une boucle supplémentaire à l'intérieur de la boucle sur les pas de temps. D'autre part, la matrice est triangularisée à chaque itération de chaque pas de temps. Le processus itératif cesse lorsqu'une certaine norme de la différence entre la solution à l'itération j et la solution à l'itération j—l est inférieure à e donné.

L'équation à résoudre est de la forme : K, (u,t)U(t )=K (u, t

1)^^

2

+ 6(10.

L a technique d'actualisation (méthode dite « du point fixe ») est coûteuse en temps de calcul, surtout pour les modèles tridimensionnels, et ce pour deux raisons : — K étant modifiée à chaque itération, il est nécessaire de refaire sa triangularisation à chaque fois : t

III -

APPLICATIONS

Les algorithmes décrits précédemment permettent de traiter les problèmes d'écoulement en milieu poreux saturé sans surface libre, et les problèmes d'écoulement en milieu poreux non saturé sans surface de suintement. Pour pouvoir prendre en compte la présence d'une surface de suintement, d'une surface libre, ou pour les calculs de nappes aquifères, les algorithmes généraux nécessitent quelques adaptations qui sont décrites dans les paragraphes III. 2 à III. 4 de cette partie.

État initial U(t ) 0

III. 1. Milieu poreux saturé sans surface libre

\ ' =0 I i = i+l ; =0

ne

j=J+l

Calcul de ^ (l/j,) Calcul du second membre K I O ^ - I + GCI-I)

Triangularisation 1

Oui

Nouveau pas de temps

Non

Nouvelle itération

Fig. 8. — Algorithme de calcul en régime transitoire (modèle n o n linéaire).

U n exemple classique est le problème de l'écoulement autour d'une palplanche. Dans ce type de calcul, les chargements se réduisent en général à des flux nuls à travers certaines surfaces; ce type de condition aux limites est satisfait implicitement sur un bord libre. Les conditions aux limites sur la charge hydraulique sont en général des conditions de charge constante sur des parties de contour, ou de charge égale à la cote (condition de pression nulle) sur des surfaces de suintement connues. Les figures 9 et 10 représentent le maillage et les réseaux d'équipotentielles et d'équipressions obtenus pour un calcul d'écoulement autour d'une galerie circulaire. L'alimentation est assurée à la partie supérieure du modèle par une condition de charge constante. L a surface de la galerie constitue la surface de suintement ou de pression nulle. Ce calcul peut-être effectué avec les programmes P P R - B I D I F ou P P R - N S A T . Dans ce dernier cas, le contour de la galerie est déclaré en surface de suintement, ce qui évite à l'utilisateur de rentrer point par point la condition h=y, cette condition étant automatiquement mise en place à la deuxième itération. L e milieu restant saturé, le processus itératif s'arrête après deux itérations. 29

Fig. 9. — Maillage pour un calcul d'écoulement autour d'une galerie.

Fig. 10. — Équipotentielles (

) et équipressions (

).

III. 2. Milieu poreux saturé avec surface libre Formulation du problème en régime permanent Considérons le cas du barrage représenté sur la figure 11. Dans le domaine Q limité par les surfaces S , S„ S , S et S , l'écoulement est régi par : h

r

h

f

- d i v ( X grâd h)=f

Fig. 11

(29)

où / correspond à un débit volumique en général nul. Les conditions aux limites sont de deux types : — celles qui s'appliquent sur des limites connues : h = h(P) sur S (charge imposée) h = h' (P) sur S' (charge imposée) ôh — k — n , = / i (P) sur S (flux imposé). dxj

sur la charge et une condition sur le flux : S, est simultanément une surface à pression nulle et une surface de courant.

h

h

iJ

!

(30)

f

— celles qui s'appliquent sur la surface libre S, et la surface de suintement S„ qui ne sont pas des données du problème : / h = z (P) sur S, et S

r

—k

iJ

—

(pression nulle)

n, = 0 sur S,

(flux nul).

L'idée générale consiste à se fixer au départ de façon arbitraire la position de S„ à lui imposer la condition de flux nul (cette condition est implicite sur un bord libre) et à ajuster sa position à chaque itération en cherchant à satisfaire la condition sur la charge h = z(P). Cela s'obtient en déplaçant les points de S, de manière à ce que leur cote soit égale à leur charge à l'itération précédente.

( 3 1 )

Toute la difficulté consiste à résoudre (29) sur un domaine qui n'est pas connu au départ, tout en satisfaisant les conditions aux limites (30) et (31). O n voit que sur la surface libre, on a à la fois une condition 30

Méthodes de résolution

Méthode de suppression d'éléments Cette méthode, qui n'est jamais employée seule, consiste à obtenir une première approximation de la surface libre en supprimant à chaque itération les éléments pour lesquels l'ensemble des sommets ont

Fig.

1 2 . - Méthode de suppression d'éléments.

Fig. 14. — Maillage initial.

- T

Fig. 13. — Ajustement final.

une charge inférieure à la cote. L a surface libre obtenue n'est pas régulière, surtout si le maillage est grossier (fig. 12). Cette première phase est en général suivie d'une phase d'ajustement consistant à déplacer les nœuds du maillage en les amenant à une cote égale à leur charge (fig. 13). Les points dont la cote est amenée à diminuer peuvent rejoindre leurs voisins immédiatement inférieurs et il faut à ce moment là supprimer les éléments « plats » ainsi générés. Cette technique est utilisée dans le groupe 06 du système R O S A L I E pour les recherches de surfaces libres en domaine bidimensionnel ou axisymétrique. Méthode de déformation d'un ensemble d'éléments Contrairement à ce qui se passe pour la méthode précédente, tous les éléments restent actifs lors du processus itératif. Le principe qui régit le déplacement des nœuds reste le même, mais on déplace non seulement les points de la surface libre, mais aussi un certain nombre de nœuds situés en dessous d'eux, en les répartissant régulièrement sur ce que l'on peut appeler des « droites mobiles ». L a limite inférieure de la zone déformable doit être choisie de façon à ne pas perturber les conditions aux limites sur la charge. Les figures 14 et 15 représentent le maillage d'un barrage avant et après déformation. Cette technique est utilisée dans le programme F L O T du système PPR. Traitement des points de ruissellement Les points de ruissellement éventuels sont situés à l'intersection de la surface libre et des surfaces de ruissellement. Comme ils appartiennent à ces dernières, on leur applique la condition h=z, ou condition de pression nulle. Comme ils appartiennent aussi à la surface libre S,, il est évident que le traitement appliqué aux autres points de S, (faire z = h) ne les fera pas bouger. Pour éviter ce blocage, on recalcule la

Fig. 15. — Maillage déformé.

Fig. 16. — Détermination de R par extrapolation.

position des points de ruissellement après déplacement des autres points, par extrapolation linéaire ou parabolique à partir des deux ou trois derniers points de S, (fig. 16). Exemple de calculs avec le programme FLO T Le programme F L O T du système PPR utilise la technique de déformation d'un ensemble d'éléments. L a figure 17 représente le résultat d'un calcul d'écoulement à travers un barrage homogène.

Limitations des méthodes précédentes Dans la mise en œuvre pratique des méthodes décrites ci-dessus, on peut se heurter à deux types de difficultés en milieu non homogène : — celles qui sont dues à une trop grande hétérogénéité des perméabilités ; — celles qui sont dues à l'obligation de respecter la géométrie des interfaces. U n trop grand contraste de perméabilité dans le cas d'un noyau peu perméable (fig. 18) obligera les éléments situés juste en aval du noyau à se déformer beaucoup, car la surface libre subit une chute rapide 31

Fig. 1 7 . — Écoulements à travers un barrage homogène (FLOT).

Cas du régime transitoire

k

k

x

Les techniques de déformation de maillage sont d'emploi délicat en régime transitoire, car il faut pouvoir traiter des mouvements importants de la surface libre, tant en descente qu'en remontée. Ces difficultés étant supposées résolues, on peut utiliser des formules de récurrence du type Galerkin (28), où l'on itère sur la géométrie à l'intérieur de chaque pas de temps; de plus, on doit à chaque pas de temps prendre en compte dans le chargement le flux qui a traversé la surface pendant ce pas de temps.

k

2

3

Fig. 18. — Utilisation d'un point de ruissellement fictif.

dans le milieu perméable. Cela a souvent pour conséquence d'entraîner des perturbations dans la partie de la surface libre située dans le noyau. O n supprime cet inconvénient en calculant la position du point R* comme celle d'un point de ruissellement, par extrapolation à partir des deux ou trois points qui le précèdent. O n utilise cette technique dans le groupe 06 de

ROSALIE

pour des rapports de perméabilité

t

t

2

I I I . 3. Milieu poreux non saturé

Bases de la méthode de calcul

—

O n a vu en I que la non-saturation du milieu poreux pouvait être traduite par une non-linéarité du tenseur ôfo de perméabilité et du coefficient du terme en — ôt (équation (7)). U n algorithme permettant la prise en compte de ces non-linéarités est décrit en II. L a donnée sous forme de fonction de v|/=h—z de la perméabilité relative k et de la teneur en eau 0 permet à chaque itération d'actualiser les coefficients du tenseur de perméabilité K et le coefficient d'emmagasinement

L'autre type de difficulté correspond surtout au cas où la surface libre traverse un interface horizontal ou sub-horizontal. L a technique de déplacement des points sur des lignes sub-verticales ne peut s'appliquer aux points R et R (fig. 19), sous peine de déformer la géométrie des interfaces. Une technique possible en cours de test consiste à appliquer la méthode classique aux points autres que R et R et à calculer la position des points R et R par extrapolation à partir des points de la surface qui les précèdent et qui les suivent. 2

à des programmes très complexes, d'utiliser une méthode de maillage » variable, la zone située au-dessus étant traitée comme une zone non

k

supérieurs à 100 environ.

t

Tout cela conduit et il est préférable fixe et de « rigidité de la surface libre saturée.

r

2

Traitement des surfaces de suintement Le long d'une surface de suintement, les conditions aux limites changent de nature au passage du point de ruissellement, l'intersection de la surface de suintement et de l'équipression correspondant à la pression nulle. Au-dessus du point de ruissellement, on a une condition de flux nul, en dessous, une condition de pression nulle se traduisant par h=z. Lors des itérations permettant l'actualisation des coefficients, les conditions aux limites sur les surfaces de suintement sont modifiées selon la méthode suivante : Première itération

Fig. 19. — Traversée de limites de zones subhorizontales.

32

O n impose la condition de flux nul sur tous les points de la surface de suintement.

Passage de l'itération i à l'itération i+1 • Cas des points à condition de flux nul à l'itération / : — si h>z, on imposera h=z; — si h_z, on maintiendra la condition de flux nul. • Cas des points à condition h=zk l'itération i : — si le flux est rentrant, on imposera la condition de flux nul. — si le flux est sortant, on maintiendra la condition h=z.

Fig. 20. — Infiltrations dans un massif.

Pour éviter certaines instabilités, on ne modifiera la condition aux limites que sur un seul point de chaque surface de suintement à chaque itération. Cette technique est employée dans le programme N S A T du système P P R .

Exemples de calcul en régime permanent Le premier exemple (fig. 20) concerne un problème d'infiltration dans un massif à partir d'une alimentation sur une partie de la surface [Vauclin et al., 1979]. Les fonctions retenues pour représenter les variations de k et de 0 en fonction de v|/ = /i—z sont de la forme : r

lsii|/ 0 =

k

'

=

P—— si xtr<0

X

(

( Y+

e si vh^o s

M

8

Dans l'exemple traité, on a pris : a = 3.10

6

p= 4 y = 4.10 5 = 2,9

4

-

1,5

-

1.

-0,5

0

Fig. 2 1 . — Variations de kr et de 8.

6 = 0,3 S

L'allure des courbes correspondantes est donnée sur la figure 21 pour diverses valeurs de p. L a courbe en pointillé représente le coefficient d'emmaÔQ

gasinement — . t3v|/ Ce type de courbe est bien adapté à des matériaux sableux ou argileux et a été introduit dans le programme N S A T du système P P R .

Fig. 22. — Vecteurs vitesses, équipotentielles et équipressions.

33

Compte tenu des importantes variations de perméabilité que l'on peut avoir au-dessus de la zone saturée, la finesse du maillage dans cette partie du modèle doit être d'autant plus grande que la perméabilité relative décroît rapidement pour des valeurs négatives de la pression. Avec des éléments de type T3 (triangles à interpolation linéaire, dont l'usage est recommandé dans le programme N S A T ) , il faut éviter que la perméabilité relative soit divisée par un facteur supérieur à 3 sur la hauteur d'un élément. L a figure 22 représente les vecteurs vitesse, les équipotentielles et les équipressions pour un flux entrant de 4 , 1 . 1 0 " m/s et une perméabilité saturée de 9 , 7 . 1 0 " m/s.

Le troisième exemple (fig. 25) concerne un barrage avec drain horizontal. Contrairement à ce qui se passerait avec le programme F L O T , le maillage ne nécessite aucune adaptation particulière.

Exemple de calcul en régime transitoire L'exemple représenté sur les figures 26 à 29 représente un calcul de barrage homogène en vidange rapide. Les valeurs des paramètres utilisés sont les suivantes : ,fc =fc,= 2.10- 'm/s x

5

5

perméabilité

a=2.10

7

P=4 Le deuxième exemple (fig. 23 et 24) reprend l'exemple du barrage homogène déjà traité avec le programme F L O T (fig. 17). Comme il s'agit d'un calcul en régime permanent, seuls les coefficients a et P interviennent dans les résultats du calcul. Comme pour l'exemple précédent, y, 8 et 8 ne servent qu'à déterminer à titre indicatif la teneur en eau en tout point.

e =o,o6 s

teneur en eau \ y = 3 980 8 =1,8

S

O n a pris a = 3 . 1 0 et P = 4,5. Une étude paramétrique a montré que pour des valeurs de P comprises entre 3 et 5, les résultats variaient peu. L a comparaison des figures 17 et 23 montrent une très bonne concordance des deux méthodes de calcul. 6

Les charges sont initialisées par un calcul en régime permanent (fig. 27). O n ramène ensuite la hauteur d'eau amont de 5 à 1 m en 8 pas de temps de 3 heures. Cette hauteur est maintenue constante jusqu'à (=10 jours. Les figures 28 et 29 représentent les résultats pour t = 30 heures et t = 10 jours.

Fig. 23. — Maillage d'un barrage homogène (NSAT).

Fig. 24. — Écoulements à travers un barrage homogène (NSAT).

Fig. 25. — Barrage avec drain horizontal (NSAT).

34

JL

Fig. 27. — Régime permanent.

Fig. 28. — Vidange rapide ( f = 3 0 h).

Fig. 29. — Vidange rapide ( f = 1 0 j ) .

III. 4. Calcul des nappes aquifères Formulation

du problème

L a détermination des écoulements dans une nappe nécessite théoriquement la résolution d'un problème tridimensionnel. Toutefois, si l'on fait l'hypothèse simplificatrice que la charge hydraulique h est constante sur une verticale (hypothèse de Dupuit), on est ramené à résoudre en domaine bidimensionnel l'équation :

div (T grad h) = C — + q ôt s

(32)

où T représente le tenseur de transmissivité. Si e (x, y) est l'épaisseur de la nappe au point M , le tenseur de transmissivité est relié au tenseur de perméabilité par la relation simple : f=e(x,y)K 35

L'épaisseur e dépendant directement de la charge h, il s'agit à nouveau, comme pour les écoulements en milieu poreux non saturé, d'un cas particulier de nonlinéarité, qui sera traité par la méthode générale décrite en IL

Technique de résolution en régime permanent Pour pouvoir déterminer l'épaisseur de l'aquifère e, il faut connaître en tout point les cotes du substratum imperméable et du toit de la formation poreuse, z et z inf

L a deuxième difficulté est d'ordre théorique : dans la formulation exposée ci-dessus, on suppose connu le tenseur de perméabilité K en tout point du modèle. Autrement dit, on considère le modèle calé au départ, et on l'utilise directement pour déterminer l'influence de diverses perturbations (infiltration, pompage, modification locale de perméabilité pour représenter un ouvrage, etc.). E n fait, les reconnaissances sur le terrain fournissent surtout des hauteurs d'eau et rarement des perméabilités.

sup

A chaque itération du processus d'actualisation de la transmissivité, l'épaisseur e est calculée de la façon suivante selon la valeur de la charge h : si / i ^ z si z

inf

s u p

, e= z


si h^z , inf

s u p

-z

(nappe captive),

i n f

e = h-z

in{

(nappe libre),

e = 0 (domaine hors d'eau).

Cet algorithme est programmé dans le groupe 03 du système R O S A L I E . L a convergence est en général obtenue pour un nombre d'itérations inférieur à 10.

Difficultés de mise en œuvre L a première difficulté est d'ordre pratique : il s'agit, une fois le maillage réalisé, de fournir pour chaque nœud du modèle la cote du substratum et l'épaisseur maximum de l'aquifère ( z —z ). sup

inf

L'utilisation d'un pré-processeur d'interpolation (programme T E R P O L ) permet de générer ces deux tableaux à partir de valeurs connues en un nombre réduit de points différents des nœuds du maillage.

Fig. 30. — Impact hydrogéologique d'une rocade (maillage).

36

L a première partie de l'étude numérique consistera donc à caler le modèle, de façon à retrouver au mieux les hauteurs d'eau mesurées avant perturbation. Ceci sera obtenu par un ajustement, souvent délicat, des perméabilités. Il s'agit de résoudre un problème du type « problème inverse », à savoir : quelle est la répartition de perméabilités qui correspond à la piézométrie observée. L a solution d'un tel problème n'est pas unique. E n effet, on peut multiplier les perméabilités par un même coefficient à l'intérieur d'un tube de courant sans changer la répartition des charges. Dans les méthodes numériques de résolution de ce problème inverse, on rajoute des contraintes d'homogénéité qui permettent de lever cette difficulté.

Exemple de calcul Il s'agit d'étudier l'influence sur une nappe de la construction d'une rocade dont certains tronçons sont construits en tranchée, entre deux parois fichées dans le substratum imperméable. Ces tronçons constituent donc un écran pratiquement étanche. Le maillage utilisé est représenté sur la figure 30. Il comporte des

Fig. 3 1 . — Calage du modèle.

CONCLUSION

Après avoir présenté un algorithme général pour la résolution de l'équation de Poisson non linéaire en régime transitoire, on a étudié quelques modifications qui permettent de traiter les cas particuliers que représentent les écoulements en milieu poreux saturé avec surface libre, les écoulements en milieu poreux non saturé et les calculs de nappes aquifères. Les techniques numériques utilisées sont celles qui ont été programmées dans les « groupes » 03 et 04 de R O S A L I E et dans les programmes F L O T et N S A T du système PPR. Les divers exemples de calculs présentés montrent l'étendue des champs d'application de ces programmes, en particulier pour les problèmes d'écoulements dans les digues ou les barrages.

'Il /

Fig. 3 2 . — Impact de l'ouvrage.

éléments quadratiques à 6 et 8 nœuds et des éléments de continuité (*) sur le tracé de l'ouvrage. Pour simuler la coupure correspondant à la partie étanche du tracé, il suffit de rendre ces éléments de continuité inactifs. Les conditions aux limites sont du type charge imposée sur A B (rivière) et B C (équipotentielle) et du type flux nul sur C D et D A . L a rivière est bordée sur une grande partie de sa longueur par une bande mince très peu perméable. L'alimentation est assurée par une infiltration sur toute la surface. L a figure 31 représente les équipotentielles du modèle calé. L a zone située à l'extrême gauche est hors d'eau. L a figure 32 représente l'effet de la construction de l'ouvrage, rendu étanche sur les parties en trait fort.

RÉFÉRENCES BILIOGRAPHIQUES GUELLEC P., ROSALIE, Système de calcul des massifs et des structures. Notice de présentation et d'utilisation, Rapp.

L C P C , 1976. ISAACS L. T., Adjustement of Phreatic Line in Seepage Analysis by Finite Element Method, Research Report N ° CE2, Université de Queensland, 1979. NEWMAN S. P., Saturated-Unsaturated Seepage by Finite Ele-

ments, J. Hydraulics Div., A S C E , 1973. VAUCLIN M . , KHANJI D . , VACHAUD G., Experimental and

Numerical Study of a Transient, Two-Dimensional, Unsaturated-Saturated Water Table Recharge Problem, Water Resources Research, vol. 15, 5,

1979.

W A N G S. Y. et al, Finite Elements in Water Resources, C. R.

3 Conf. Inter, sur l'utilisation des elements finis en hydrogéologie, Université du Mississipi, 1980. e

ZIENKIEWICZ O. C ,

The Finite Element Method in Enginee-

ring Science, Me Graw-Hill, Londres, 1971.

37

The automatic design calculation of earth dams and dykes J.L. BANGRATZ - H. JOSSEAUME - B. LEBRUN

The purpose of the DIGUE conversational programme is to facilitate the application of the NSAT hydraulic design calculation programme and the PETAL programme for calculating slope stabilization (modified under the name STADIG) in several configurations of heterogeneous dykes built on two-layer soil systems. The dyke is in four parts (upstream slope, centre, downstream slope, downstream drainage), and there are several variants of each part. Fitting these variants together makes it possible to represent quite a large number of different configurations which may have upstream and downstream berms, a core with or without vertical chimney drains, a downstream drainage blanket, a toe fill, an upstream or central cut-off wall, and downstream relief wells. Modifications to the PETAL programme relate to the definition of the potential slipping circles and the improvement of performance characteristics in seeking values of the pore pressure along the potential slipping circles. Bull, liaison Labo. P. et Ch.. n°130,

mars-avr. ¡984,

p. 5-14.

Interpretation of measurements of pore pressure in the Maurepas-Courance dyke. Anisotropy of permeability of the structure H. JOSSEAUME - C. MAHDAVI - J. JARDIN

The authors give an account of the study of the hydraulic behaviour in service of a homogeneous impounding dam 15 metres high, constructed in compacted clay of low plasticity. This study, whose main purpose was to characterize the anisotropy of permeability of the structure, was based on: - Measurements of pore pressure in the structure, using piezometric instruments set in place during its construction, and additional piezometric instruments specially installed. - Measurements of permeability in situ, made in the dyke and its foundation. - Theoretical studies of the permanent flow through the structure, using the method of. finite elements and allowing for various conditions of anisotropy of permeability. Analysis of the experimental results and comparison of the theoretical and experimental values of pore pressure reveal an almost isotropic behaviour of the dyke in respect of permeability. Bull, liaison Labo. P. et Ch., n°13l,

Examples of specific problems relating to low dykes C. M1EUSSENS

mars-avr. 1984, p. 15-24.

Observations on the Maurepas-Courance dyke A. DELFAUT - J. JARDIN - R. BALDIT

The Maurepas-Courance storage basin is part of the development programme of the new town of Saint-Quentin-en-Yvelines. It is a storm water reservoir serving as a buffer basin for the rainwater of the plateau, whose urbanisation has increased the peak flow of flood water at the bottom of the valley. The task assigned to the LROP was more particularly to keep track of the behaviour of the body of the dyke during its construction, and subsequently after it had been put into service. For this purpose it was necessary to fit two transverse profiles of the dyke with pore pressure cells, settlement recorders and piezometric tubes. Optical measurements were made on concreted reference points so as to assess movements of the body of the dyke. The authors report the principal results obtained during the construction of the dyke, during pondage, and during the operation of the basin. To attempt to understand the real hydraulic functioning of the structure, it appeared worth-while to compare the observations with the initial theoretical analyses carried out in the context of research on the anisotropy of the silts used in the body of the dyke. Clear conclusions were not reached, and it is difficult to take account of the complexity of real phenomena in theoretical studies, the more so in the planning stage. Bull, liaison Labo. P. et Ch., n°131,

Observations on the Roussières dyke J. JARDIN - R. BALDIT - A. DELFAUT

The author relates some examples of specific problems encountered during the design or the follow-up after putting into service of five dykes located in the area covered by the Laboratoire Régional des Ponts et Chaussées in Toulouse. 1. The case of a dyke built on excessively permeable soil. The dyke and underlying soil waterprofing was realized by means of a grout curtain which proved inadequate in view of the low admissible leakage flow. 2. A low dam holding back a very considerable quantity of water, whose downstream face was shown by appropriate instrumentation to be tending towards certain failure. 3. Investigative methods employed to study the permeability of the soil underlying a dam, which was a priori geologically very unfavourable because of the hazardous presence of permeable horizons. For the body of the dyke, the influence of lime treatment on the permeability of the compacted material was measured. 4. The design of flood control dam so that it could best ensure maximum storage. The problems of internal erosion and by-passing are more particularly examined. 5. The special case of an expanse of water located near a river, necessitating a study of the stability of the sealing curtain should the dyke function in reverse (pond empty, river in flood). This same example also includes a description of a simulation of plane flow by electrical analogy. In conclusion, in giving an account of sometimes simplistic methods of study, the author raises the problem of adequation between the level of studies and the economy of the project. Bull, liaison Labo. P. et Ch., n°130,

mai-juin 1984, p. 23-28.

mai-juin 1984, p. 5-22.

Located on the territory of the new town of Saint-Quentin-en-Yvelines (in the Yvelines département) the Sourrières dyke is a dam of medium height (about 15 metres) holding back a 400,000 m rainwater reservoir. It is a homogeneous silt dyke with a drainage system (drain shaft and drainage blanket). Though of relatively modest dimensions, the structure, which concerns public safety, was equipped with important inspection facilities. Observations of its hydraulic and mechanical behaviour were made duringfillingof the reservoir to the highest level on completion of its construction, and for two years at the permanent level after the dyke had been put into service. In particular, after additional piezometers had been installed a study was made, at the permanent level, of the hydraulic anisotropy of the silt constituting the body of the dyke. After reviewing the studies carried out at the planning stage, the authors report the results of these observations. The almost isotropic behaviour (k /kv = 1) of the silt is the most important result of these observations. The authors believe that this fact is mainly linked with the way the material was applied. 8

n

Bull, liaison Labo. P. et Ch., n°131,

mai-juin 1984, p. 29-44.

The construction and behaviour of an earth dyke at high altitude M. LAGRANGE - J. PERR1N - J. JOSSEAUME

The provision of a storage reservoir to improve the water supply of a winter sports resort necessitated the construction of a low dyke (9.50 m) at high altitude. After presenting the principal geotechnic characteristics of the site, which led to the construction of a homogeneous dyke, and a description of the dyke itself, the authors give an account of the problems which arose during and after the work; problems connected with the laying of a geotextile between the body of the dyke and the protective enrockment, and considerable leakages when the reservoir was filled, necessitating the injection of the rock foundation. The article also reports measurements of pore pressures in the dyke, and their interpretation. Bull, liaison Labo. P. et Ch., n° 132, juil.-août 1984, p. 13.

Application of the method offiniteelements toflowin a porous medium P. HUMBERT

After a reminder of the basic equations governing flow in a porous medium, a technique of solution, using the method of finite elements, of Poisson's non-linear equation under transitory conditions is presented. The adaptations to be made to this general algorithm in order to solve certain special problems are explained and illustrated by examples in the case offlowin unsaturated porous media and in the case of water-table calculations. Bull, liaison Labo. P. et Ch., n" 132, juil.-août 1984, p. 21.