UNCONFINED FLOW T H R O U G H POROUS MEDIA C l a u d i o Baiocchi, V a l e r i a n o C o m i n c i o l i , 1 and U g o M a i o n e 2

S O M M A R I O : Viene esposto ten nuovo metodo di studio di problemi di fi#razione a superficie libera basato sull'utilizza~one di disequagioni variazionali the permette di ottenere teoremi di esistenza e di unicit~. Vengono introdotti nuovi procedimenti di calcolo numerico che forniscono la soluzione del problema come limite di una successione monotona cres:ente e "di u#a monolona decrescente. Tali procedimen*i, oltre che essere giustificati dal punto di vista della convergenza e della stabilita, risultano competitivi con quelli attualmente in uso, sia per la rapidith di convergenza, che per la semplicita di programmazione. Vengono riporlati alcuni risultati numerici. S U M M A R Y : A new method of studying unconfinedflow through porous media has beet, illustrated. This method is based on the use of variational inequalities which allows us to obtain uniqueness and existence lheorems. New numerical methods have been adopted which give the solution of these problems as limit of an increasing and decreasing monotonic sequence. Such procedures, other than being justified from the point of view of convergence and stability, are competitive with those actually in use both for the rapidity of convergence and for the simplidty of the programmation. Below, we shall illnstrate some of the numerical results.

in which:

u=y+

P-Y

denotes the "piezometric head." This new theoretical method has allowed us to obtain existence and uniqueness theorems and besides has originated new numerical methods. The main features of such methods are the rigorous theoretical base on which they are founded and besides the simplicity of programming and the speed of execution for which they compete very well with the methods already known. Before examining in detail these problems, it is useful to recall some of the basic notions concerning the flow of fluids through porous media. Let us consider the classic case of the flow of an incompressible liquid through an earth dam which separates two reservoirs of levelsyl and y~; the floor of the dam is supposed impervious (fig. 1).

y y = Y ( xl (x)-

E r~,

1. Introduction. For some years, in collaboration with the Istituto di Idraulica della Universit~ di Pavia, the Laboratorio di Analisi Numerica del Consiglio Nazionale delle Ricerche (CNR) has conducted a theoretical study of unconfined flow through porous media which has been largely referred to in specialized scientific papers (see, for example, [1], [2], [4], [5], [6]). The most recent results of such research are the object of this article. The investigation was originated by Baiocchi [1], who has given to such movement processes a new mathematical formulation in terms of variational or quasi-variational inequalities, by introducing the new unknown function:

z('¢,:) =

0

c

Denoting by D the dam and by at2 the flow domain, from the mathematical point of view the problem is described by the equations:

0~

Au=O

in.q A = ~ + ~

02 )

u l =0 .It, = :

&

(1.1)

(1.2) 0.3)

0.4)

u Ir, = Y

SEPTEMBER 1975

x.

Fig. 1

t - - u(,¢, t)] dt

z Istituto di Matematica, Universit~ di Pavia e Laboratorio di Anatisi Numerica del CNR. 2 Istituto di Idraulica, Universith di Pavia.

b

0

0

denotes the outward normal deriv.)

0.5) 151

where I ~, indicates the "seepage face" and I', the free line. We recall also that (1.1) follows from the assumed validity of Darcy's law:

and only ooe solution Z. Thus it makes sense to consider the fuoctional transformation from W to the corresponding Z- Such transformation will be denoted as ~J" (i.e., Z = ~-'W). We will also introduce the transformation ,9° defined by:

V = - - K grad u

(Sail;") (x) = ( J ' W ) (x, Y(x)), a <~ x ~ b

where V is the seepage velocity and K is the coefficient of permeability. The other relations define the following boundary conditions:

Returning to problem (1.1) . . . . , (1.6) let us introduce a change 'l of unknowns setting:

a) the piezometric head is constant along A F and BC; b) the pressure is zero along the free line Ft and the seepage face 1;; c) the seepage velocity is directed tangentially to the base (which is impervious) and to the free line (which is a streamline). For the developments which follow, it is of interest to note that since the fluid emerges from the porous medium at the seepage face I ; , we have the following condition: 0u --~Tni.% ~< 0

Z(x'Y) = .i~ [t--/~(x, t)] dr, for (x,`),) ~ D--,

(2.3)

(2.4)

where t7 is the extension of u which equals), outside .(2. Setting: g ( x ) = ~:(x, l'(x)) for a ~ x ~
(2.5)

(1.5) it follows:

W is a fixed point of Sa, that is W - - S a i l ;',

(2.6)

0.6) while (1.6) becomes:

From a mathematical point of view, the problem now described is a classical free boundaO, problem and consists in determining a curve y = q,(x) (the graph of which is 1~ U l'l) joining F to C; in the "subgraph" /'2 of q, we must solve a boundary value problem with "too many" conditions ((1.4) and (1.5)) on the "free part" I'z of OO. For a mote extensive description of the seepage problems see, e.g., [7] and the Bibliography there quoted; for some recent developments see, e.g., [11].

(W(x))"<~ 0 on ]o, c[

Vice versa, knowing a function W which satisfies (2.6), (2.7) we can deduce a solution of the free boundary problem by setting: o = {(=,y) e D 1(.~,-~3(x,y) > iv(=)}; (2.8)

,,(x,0,) = y 2.

(2.7)

( . ~ g ) (xa,) for (x,O') ~ ~2

M a t h e m a t i c a l treatment.

In order to study problem (1.1l . . . . , (1.6) it is useful to first examine the following auxiliary problem. Let W(x) be any function defined on [a, hi, then we look for a function Z defined o n / i f ( = closure of D), which satisfies the relations

0 <~ A z ~ 1 o n D ;

ZIz = O where z(x,y) < W(x);

A z = 1 where Z(x,`),) > W(x) a

Remark 2.1. The idea of solving a free boundary problem by studying the fixed points of a suitable functional transforrnatioo is already present in many papers, especially of numerical type, devoted to the subject. For example, let us consider the transformation defined as follows: at each "smooth" curve ~(x) joining F to C and with 9~(x) ~< ~< Y(x) we associate the corresponding subgraph ~w; denoting by l'z and l', respectively the curves {(x, ~(x)l

19'(") < Y(=)} and {(,,, V,(~M(~)= Y(=)}, tet ,,,~ be (2.1)

z l ~ = o; ~1,'~ =`)'-`),~; z~l~ = o; ~,,Ic.'~ =`),-`),2 (2.2)

the solution of problem defined by (1.1), (1.2), (1.3), (1.5) with Q = / 2 o. Setting ~ = u~(x, ft.(x)), (1.4) will be verified if and only if q~ is a fixed point of Y. However, the study of the transformation ~, and of its fixed points is not easy because ~'~ is obtained by means of the solution

Under wide assumptions of regularity on the curve y = Y(x) problem (2.1), (2.2) has, for each fixed W, one

3 These relations are usually summarized as follows: A z ~ H(Z--IF') where H is the multivalued function defined by H(I)= {0} for I < 0; H(1) = {1} for I > 0; H(0) = {hl0 ~
152

,t Introduced in [2] and which generalizes the one introduced in [4]. Such transformation has been found useful in various free boundary problems. Besides to the seepage Froblems, it has been applied successfully to fluidodynamic problems and problems of the Stefan type; cf. [3] and the related Bibliography. 5 Proofs will be given in a forthcoming paper.

MECCANICA

of a boundary value problem in a domain which depends i/self On the contrary, the study of the fixed points of 5 a involves W in the equation, but with respect to the fixed domain D, which greatly simplifies the treatment.

for each n, IV, is a fixed point ofoq'; the sequence IV,+ converges increasing to a function W®, which is a fixed point of.9 °. In addition W® satisfies (2.7) (2.14) and every other W fixed point that satisfies (2.7) is such that IV>~ W®.

Remark 2.2. Problem (2.6), (2.7) can be considered as a

Thus the function W~o gives the minimum solution of the problem.

on the ~o.G

quasi-variational problem (cf. [2]).7 It will, in particular, follow the existence of a "maximum solution" and of a "minimum solution." Following the technique introduced in [8], we can prove that the iterative procedure defined by: W°(x) = 0; W,,+~ = S,"tU, (n -----0, 1, 2 . . . . )

(2.9)

Remark 2.8. Naturally, we have W e ( x ) ~ W°~(x). The problem of the coincidence between W,o and IV= (and thus of the uniqueness of the solution of (1.1) . . . . (I.6)) remains, in the general case, an open problem. However, the property can be proved, for instance, if it is assumed that D has "the right wall vertical." Under such a hypothesis (2.7) can be specified as follows: ( I V ( x ) ) " = 0, 2

is such that: W, converges decreasing to a function ll7=, which is a fixed point of ~9v and every other eventual (2.10) IV fixed point of &v is such that IV ~< IV%

that is W ( x ) = q ( x - - a ) - - y z / 2 with q = discharge. In the general case the numerical experiments carried out (cf. section 3) suggest the conjecture that in any case we must have : IF/'~ = IV=. 3.

In addition:

N u m e r i c a l results.

(2.11)

The procedures now described, besides proving the existence of IV~ and W® are of a constructive 0,pc: to any discretization of problem (2.1), (2.2) (that is, of the functional transformations .Y, ~ we can associate the sequence

and thus It/°: is the maximum solution of (2.6), (2.7). In order to obtain the minimum solution we shall define, firstly, Zo as the solution of the problem (analogous to (2.1), (2.2), but simpler because it is linear):

IV~ (h = discretization parameter) analogous to the one introduced in (2.9). Also, with a discretization of the operator "cone" (of. Remark 2.2) we can define a sequence Wn.,. Formulas like (2.8) will then give the discrete sets

(IV%~))" ~< o on

]0, ~[

II 0 = 1 on D ; boundary conditions (2.2).

(2.12)

We can prove that: IV0(x')--Z0(x , Y(x)) is a fixed point ofE~° and is less than or equal to every other fixed point; however, (2.7) is not in general satisfied by 1V0. Let us introduce now a procedure of "concavization" on ]0, el: for any function f defined on [a, b] we shall indicate with cone ( f ) the function that on [a, 0] and on [c, b] coincides with f ; while on ]o, c[ it is concave and is less than or equal to each concave function greater t h a n f , s With this notation we put: w0(x) = ~0(x, Y(x)); g,,+~ = 5~ (eonc (II<,)) (2.13) (n = 0, 1, 2 , . . . ) Now we can prove that:

s From the numerical point of view a sequence of boundary value problems must be solved in a family of domains each one depending on the solution of the preceeding problem. Quasi-variational inequalities were introduced in [8] for the study of impulse control ; also in this field they describe free boundary problems. a On ]o, el, cone (f) solves the problem of the "elastic line forced under an obstacle." From the mathematical point of view it corresponds to a simple minimum problem. SEPTEMBER 1975

/'22 and On,,, and the approximations of the function u and of the free boundary. Bearing in mind that problem (2.1), (2.2) is equivalent to a variational inequality, the specialized literature (cf., for instance, [9], [12], [13], and the References of these papers) gives for such a problem various types of discretizations (finite differences, finite e l e m e n t s . . . ) . These are, at the same time, efficient from the point of view of the calculating speed, of programming simplicity, and in addition they provide the convergence of Wnn to IV", of Z~ tO Z, etc. Referring to [10] and to works in preparation for a more detailed analysis of the discrete problem and of the results of the numerical experiments carried out, we give a brief description of the algorithms used in order to obtain lVt~ and lVn,~. Making use for simplicity of the finite differences method, we introduce in the plane (x,y) a net with the step b = = (hi, &), and we denote by Dh the set which discretizes D and by/In the classic five points formula which approximates the operator A. For each fixed b a sequence W~ is given analogously to what has been done for the continuous problem (2.9). Thus, it is necessary to introduce a discretization of problem (2.1), (2.2). In practice, for each fixed n, we must solve problems of the following type, where we indicate, in a general way, with Zh the unknown function and with ph an assigned function:

Anzn ~ H(Zn --ph) + boundary conditions that discretize (2.2)

(3.1) 153

Such a problem is a "complementarity system" and is equivalent to a variational inequality in a finite dimensional space (in some cases to a m i n i m u m problem of a quadratic functional, cf., for example, [1]). I n matrix form, it can he written:

f ( M ) ~ H(Zn(M ) - - p h ( M ) ) + A z n ( M )

(3.2)

where M ~ Dn, A is the "stiffness" matrix related to An a n d f collects the terms coming £rom the boundary conditions in (3.1). Thus problem (3.1) is a O,stem of nonlinear equations. To solve it we can use iterative methods. Since the matrix A is an M-matrix and H is a monotone operator, there is convergence of the following procedure which is a generalization of the relaxation method (SOR) already k n o w n to solve the linear systems. For brevity, we shall indicate with Zt the components of vector zn(M), i.e., i = 1, 2 . . . . . Nh ( = the n u m b e r of points of Dn) : ..rio) } we give a sestarting from an arbitrary vector {x; t (m)~ quence of vectors lZ~ j, solving the equations: f-1

(ttl)

,of, + (1

- - ,o) auz,

(m)

- - o J ~ auzj

(/n+l)

Maximum solution.

(re+l)

~a,z~

(re+l)

+ oH(z,

--PO

(3.3)

j-t+l

where ~ is a real parameter, which is suitably chosen in order to accelerate the convergence of the procedure. For each fixed i it is necessary to solve again a nonlinear equation. However, from the definition of H , the corresponding solution can be easily given in an explicit form. I n practice, this involves the search of the m i n i m u m of the following function:

1 t-p, t - + j ( / ) = - ~ - aut 2 + ~o (m)

+(I--cO) Zl

t-1

--0)

+ !t-p,I 2

(w fi +

(m+l)

~ a,SZ j

--co

~1 a u z j

1) one iteration of the SOR method for each fixed n; outer iterations: 158; execution time: 35 sec.

0 = maxlW~÷l(ib0 -

IV~(ih,)l = 0.205 10 -5 ( f o r

n = 157)

(m +1)

= max [Zn

(ihl, _]h,,). - - Z~(")_(ih1,fl.~.)l __; = 0.38 1G-'

¢,t

discharge q = 0.31503 2) for each fixed n a maximum number of 10 iterations of SOR method. outer iterations: 137; execution time: 107 sec.

--

J-1 Nh

for each fixed h, the sequence IV~ is monotonically decreasing, while lVn.n is monotonically increasing. This permits us to numerically enclose the solution. Another aspect of the procedure, which seems useful to point out, is the fact that in the outer iterations for calculating IV~ and lVh~ the discreatization domain does not change. This last feature is important from the programming point of view. In order to show the efficiency of the method we shall now report some numerical results for a particular problem. This calculation was carried out in double precision, on the Honeywell 6030 of the Centro di Calcolo Numerico della Universit?t di Pavia. The problem refers to a dam with a trapezoidal section defined by the points: (0,0), (2.5, 0), (1.5, 1), (1, 1). The two levels a r e y i = 1,3,.o ----0.2 Discretization (finite differences) : hi ----h2 = 1120 Relaxation parameter (in (3.3)): co = 1.85.

- - Co ~ a,lz s 1=1

Nn

Keeping in mind the significance of the concavization operation and above all the dimension of vector Wn,, (compared to the one of vector Znm) we can use some simple direct algorithms (for example, we successively raise the nonconcave "parts" by straight lines). As in the continuous problem also in the discrete one,

0 = 0.7 10-5 (for n = 136); e = 0.37 10-4; q = 0.31955

(m). ) t.

J=l+l

Minimum solution. Then, for the calculation of Z~, we must carry out an outer iteration (in n) and an inner iteration (SOR). In practice, we can suitably combine the two iterations: for example, one iteration of SOR for each fixed n. Such a variant keeps (see following numerical results) the convergence and permits, on the other hand, a significant saving in execution time. The cost, in this way, for the

The SOR method is stopped for each fixed n when:

evaluation of Zff, is equivalent to solving one variational inequality. I n order to calculate Zn~, there is a procedure analogous to the one used for Z~- It is only necessary to add, for each fixed n (of. (2.13)), the concave envelope of lVn.n.

m a x [ c o n c (IP'h..(ih,)) - - ~ , , , , ( i h l ) ] = 0.28 10 -4 ( f o r n = 14)

] 54

e ~< 0.4 10 -s outer iterations: 14;

execution time: 125 sec.

t

q = 0.31326. The following result is interesting with respect to the MECCANICA

theoretical conjecture that It7¢o= IV® (c£ Remark 9.3). r h e difference in the maximum n o r m between IVff (calctflated as in 2)) and IVn,, is = 0.87 10 -2. Making hi = =- h.o -----1/30, and keeping the same stopping tests for the SOR method and for the search of the fixed point, such a difference becomes: = 0.56 10-L The following table shows the flow net related to the m i n i m u m solution (coinciding, practically, with the one related to the maximum solution as the values, practically coinciding, of the corresponding discharges shown). 2nd National Congress, October 1974

REFERENCES 1 BAIOCCHI, C. Ann. Mat. Pura e AppL, (IV) XCII (1972), 107-127. Nota preliminare su C. R. Acad. Sc., Paris, t. 273 (1971), 1215-1217. 2 BAIOCCHI, C. C. R. Aead. Se. Paris, t. 278 (29 Avril 1974) 1201-1204. 3 BAIOCCHX, C. Corso C.I.M.E. II Ciclo, Bressanone 1973, Edizioni Cremonese, Roma, 1974. 4 BAmccHL C., COM~NCmLI, V., MAGENES, E., PozzI, G. A., Ann. Mat. P,ra e AppL, (IV) XCVII (1973), 1-82. 5 BAIOCCHI, C., COMINCIOLI, V., GUERRI, L., VOLP[, G., Caleolo, X Fasc. 1 (1973), 1-86. 6 Bxmccru, C., MAtONE, U., XIII Convegno di Idraulica e Costruzioni Idrauliche, 21-23 sett. 1972, Milano. SEPTEMBER 1975

FLON NET

Fig. 2

7 BEAa, J., Dynamics of fluids in porous media, Els. Publ. New York, 1972. 8 BENSOUSSAN,A., LIONS, J. L. C. R. Acad. Sc. Paris, t. 276 (1973), 1189-1193, 1333-1337. 9 CEA, J., GLOWINSKI,R., Int. J. Comp. ll4ech., vol. B (1973) 225-255. 10 COMINCIOLI,V. Pubblicazione n. 79 del L.A.N, del CNR, Pa~i,, 1974. 11 DATEI, C., GrmTTI, A. Conferenza al XIII Convegno di Idraulica e Costruzioni Idrauliche, 21-23 sett. 1972, t~[ilano. 12 GLOWXNSKq R., LIONS, J. L., TREMOLmRES, R., Risolution numlrique des in~quations de la l~lieanique e de la Phytique. To be published by Dunod, Parigi. 13 Mosco, U. Corso C.I.M.E. II Ci¢to Erice 1971, Edizioni Cremonese, Roma, 1973. 15"5