Homogenization of the layer-structured dam problem with isotropic permeability

Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

Contents lists available at SciVerse ScienceDirect

Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa

Homogenization of the layer-structured dam problem with isotropic permeability Sébastien Martin a,∗ , Carlos Vázquez b a

Université Paris-Sud, Laboratoire de Mathématiques, CNRS-UMR 8628, Faculté des sciences d’Orsay, Bâtiment 425, 91405 Orsay cedex, France b Universidade A Coruña, Dep. Matemáticas, Facultad de Informática, Campus Elviña, 15071-A Coruña, Spain

article

info

Article history: Received 10 September 2012 Accepted 9 April 2013

abstract In this paper we address the homogenization of a free boundary model for a stratified porous dam problem. The motivation of the problem lies in the presence of a highly oscillating permeability to take into account the multilayer structure. Although several papers in the literature have partially studied this problem, we provide new additional results by using the periodic unfolding method. Furthermore, by using appropriate numerical techniques based on the combination of characteristics, duality and finite elements methods, several examples illustrate the here stated theoretical results and show the performance of the numerical methods. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction This paper deals with the homogenization of the heterogeneous dam problem, as formulated by Brezis, Kinderlehrer and Stampacchia [1] and Alt [2]. This study is motivated by the varying permeability which allows to take into account a periodic layered structure of the porous medium and, as usually in the homogenization technique, we aim at describing the resulting homogeneous dam problem. This has been partially studied by Rodrigues [3] and, here, we provide additional results by using the periodic unfolding method [4]. Numerical computations based on the extension of the techniques proposed by Bermúdez and Durany in [5] to the case of variable and nondiagonal permeability matrix allow to illustrate the obtained theoretical results. The dam problem has first been stated using variational inequalities by Baiocchi and coauthors [6–8] and Benci [9], although this approach is only suitable for dams with vertical walls (typically rectangular dams). The formulation of the dam problem for domains with general shapes has been introduced by Brezis, Kinderlehrer and Stampacchia [1] and Alt [2]. In order to recall this formulation, let us consider a bounded domain, denoted Ω , with a locally Lipschitz continuous boundary ∂ Ω = Γper ∪ Γ0 ∪ Γa (see Fig. 1). Here, Γper denotes an impervious boundary, Γ0 is the boundary in contact with the open air and Γa is the boundary in contact with water. The basic problem is to find the pressure p and the fluid saturation θ in Ω . Introducing the permeability of the porous medium, denoted by k, and the gravity g := −g e2 , the strong formulation, based on Darcy’s law (see [10] for historical references), is given by the following set of equations: div (θ k g − k∇ p) = 0,

∗

p ≥ 0, θ ∈ H (p),

Corresponding author. E-mail addresses: [email protected] (S. Martin), [email protected] (C. Vázquez).

1468-1218/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.nonrwa.2013.04.002

(1)

2134

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

Fig. 1. Dam domain.

where H denotes the Heaviside graph. The following boundary conditions are considered:

 p = α,   p = 0,  (θ k g − k∇ p) · n ≥ 0, (θ k g − k∇ p) · n = 0,

on Γa , on Γ0 , on Γ0 , on Γ ,

(2)

with, for example,

α(x) = hi − x2 ,

on Γa,i ,

where hi denotes the water level in the ith reservoir. Here, n denotes the unit normal vector that point outward from Ω on the boundary ∂ Ω . Remark 1. In the strong formulation, one should notice the boundary condition on Γ0 corresponding to a sign condition on the fluid flow: it is designed to eliminate the non physical solutions and means that no water flows into the dam through Γ0 . The formulation of the dam problem for general shapes was introduced by Brezis, Kinderlehrer and Stampacchia [1] and Alt [2]. The variational formulation of that problem is given by:

 ∞ Find (p, θ ) ∈ V0 ∩   Wα × L (Ω ) such that:  k ∇ p · ∇φ − θ k g · ∇φ ≤ 0, ∀φ ∈ V+ ∩ W0 , (A)  Ω  Ω p ≥ 0, θ ∈ H (p), a.e. in Ω , with the following definition of the functional spaces:





V0 = v ∈ H 1 (Ω ), v|Γ0 = 0 ,





Wα = v ∈ H 1 (Ω ), v|Γa = α ,





V+ = v ∈ H 1 (Ω ), v|Γ0 ≥ 0 ,





W0 = v ∈ H 1 (Ω ), v|Γa = 0 .

In the case of a homogeneous porous medium (k ≡ 1), existence was proved after passing to the limit on a penalty parameter η, as the saturation function θ was approached by a regularized nonlinearity Hη (p) that approximates the behavior of the Heaviside graph. Uniqueness of the so-called Γa -connected solution was established by Carrillo and Chipot [11] and Alt and Gilardi [12]: roughly speaking, if the shape of Ω is such that no pool appears, then problem (A) admits a unique solution. Several authors also addressed the regularity of the free boundary, as function θ is a characteristic function, meaning that the free boundary divides full-saturated regions ([p > 0] or [θ = 1]) and full-dry regions ([p = 0] or [θ = 0]). In particular, the regularity of the free boundary ∂[p > 0] was investigated by Alt in [13] proving that it is an analytic curve y = Φ (x) when Ω is a Lipschitz domain. In the case of a heterogeneous porous medium, the formulation has already been described by Alt in [2]. The penalty method remains valid in order to prove existence of a solution. Stavre and Vernescu [14] and Friedman and Huang [15] extended previous results described for the homogeneous case in the following sense: when k is weakly nondecreasing in x2 , there exists a unique Γa -connected solution, the saturation function is still a characteristic function, and the free boundary is continuous. The numerical computation of solutions to the dam problem can be performed by several methods: Alt [16] and Marini and Pietra [17] built a numerical method based upon a discrete analogue of the continuous problem, by introducing finitedimensional spaces and a fixed-point procedure. Bermúdez and Durany [18,5] proposed to solve a transient version of the dam problem, using a combination of the method of characteristics and the finite element method; then the solution of

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

2135

the nonlinear discretized problem is obtained by using a duality iterative algorithm. Note that Chakib and coauthors [19] also proposed to tackle a transient version of the dam problem: using a time discretization, the steady-state solution is considered as an optimal shape design formulation on a transient version of the dam. Recently, Nazemi and coauthors [20] proposed a numerical technique based on an optimal control approach of the problem, governed by variational forms on a fixed domain; then by using an embedding method, the class of admissible shapes is replaced by a class of positive Radon measures and the optimization problem in measure space is then approximated by a linear programming problem. However, their method only applies when the saturation function is a characteristic function. For the sake of clarity, let us introduce some related formulations of the dam problem, based on different model assumptions:

• Penalized dam problem. Let η > 0, and let us define  0, if p ≤ 0, Hη (p) = p/η, if 0 < p < η, 1, if η ≤ p. Obviously, when η goes to zero, p → Hη (p) tends to mimic the behavior of the Heaviside graph p → H (p). Thus we can introduce the associated penalized problem:

 Find pη ∈ V0 ∩ Wα such that:   k ∇ pη · ∇φ − Hη (pη )k g · ∇φ ≤ 0, (Aη )  Ω  ηΩ p ≥ 0, a.e. in Ω .

∀φ ∈ V+ ∩ W0 ,

Passing to the limit with respect to η in the penalized dam problem allows to get the existence result for the dam problem. • Dam problem with leaky boundary conditions. From a physical point of view, one may have to prescribe the flux on Γa , instead of the value of the pressure. More precisely, the flux may be governed by a function of the jump of pressure across the boundary. This kind of conditions are often referred as leaky boundary condition, see [21], and the Dirichlet boundary condition is then replaced by

(θ k g − k∇ p) · n = −β(·, α − p),

on Γa ,

where β is assumed to satisfy some regularity assumptions: (i) x → β(x, 0) ∈ L2 (Γa ); (ii) x → β(x, p) is measurable for every p ∈ R; (iii) ∃C > 0, |β(x, p1 ) − β(x, p2 )| ≤ C |p1 − p2 |, for almost every x ∈ Γa , for all (p1 , p2 ) ∈ R2 ; (iv) p → β(x, p) is nondecreasing for a.e. x ∈ Γa ; (v) β(x, p) ≥ 0 for a.e. x ∈ Γa , for all p ∈ R. The related variational formulation is then:

 ∞ Find (p, θ ) ∈ V0 ×   L (Ω ) such that:   k ∇ p · ∇φ − θ k g · ∇φ ≤ β(·, α − p) φ, (B)  Ω Γa  Ω p ≥ 0, θ ∈ H (p), a.e. in Ω .

∀φ ∈ V+ ,

In the homogeneous case, Carrillo and Chipot [22] proved the existence and uniqueness of the Γa -connected solution. The results were extended to the heterogeneous case by Chipot and Lyaghfouri [23] and Challal and Lyaghfouri [24–26]. The rest of the paper is organized as follows. Section 2 contains the main theoretical results concerning the homogenization of the problem. Section 3 describes the numerical methods to approximate the solution of the small parameter dependent problem and the homogenized one. Finally, some appropriate numerical examples to illustrate the theoretical results and to validate the performance of the numerical methods are presented in Section 4.

2. Homogenization of the problem Let us assume that the behavior of the porous medium is stratified. Such a multilayer structure of the medium can be modeled by a highly oscillating permeability, which can be described by:



kε (x) = K x,

x

ε

,

∀x ∈ Ω .

We suppose that kε is uniformly bounded and coercive (with respect to ε ). We denote by (Aε ) the related dam problem with oscillating coefficients

2136

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

 ∞ Find (pε , θε ) ∈ V0 ∩   Wα × L (Ω ) such that:  kε ∇ pε · ∇φ − θε kε g · ∇φ ≤ 0, ∀φ ∈ V+ ∩ W0 , (Aε )  Ω  Ω pε ≥ 0, θε ∈ H (pε ), a.e. in Ω . We aim at describing the behavior of (pε , θε ), as ε goes to 0, and characterize the limit functions as a solution of a so-called homogenized dam problem whose mathematical structure is close to the one of the initial problem. Although the core of the asymptotic analysis will not focus on the formulation including leaky boundary conditions, the homogenization method can also be applied and we will just provide the homogenization result in the case of leaky boundary conditions, underlining the main differences with the case of Dirichlet boundary conditions. 2.1. Preliminaries: the periodic unfolding method First we recall some useful definitions and results for the periodic unfolding method which has been introduced Cioranescu and coauthors [4,27]. This method has strong links with the two-scale convergence technique which was introduced by Nguetseng [28], and further developed by Allaire [29], Lukkassen and coauthors [30]. Lemma 1. The separable Banach space L2 (Ω ; Cper (Y )) is dense in L2 (Ω × Y ). Moreover, if f ∈ L2 (Ω ; Cper (Y )), then x → σε (f )(x) = f (x, x/ε) is a measurable function such that

    σε (f ) 2

L (Ω )

    ≤ f  2

L (Ω ;Cper (Y ))

.

Let us now briefly introduce periodic unfolding methods [4]: Definition 1 (Unfolding Operator). For any x ∈ R2 , let be [x]Y the unique integer combination

n

j =1

kj bj of the periods such

that x − [x]Y belongs to Y . We also define {x}Y = x − [x]Y ∈ Y , so that, for each x ∈ R2 , one has

  

  x=ε

x

ε

+ Y

x

ε

. Y

Then the unfolding operator Tε : L2 (Ω ) → L2 (Ω × Y ) is defined as follows: for w ∈ L2 (Ω ), extended by zero outside Ω ,

    x + εy , ε Y

Tε (w)(x, y) = w ε

for x ∈ Ω and y ∈ Y .

Let us now recall main results by Cioranescu, Damlamian and Griso [4]: Proposition 2. One has the following integration formula

 Ω

w=

 Ω ×Y

Tε (w),

∀w ∈ L1 (Ω ).

Proposition 3. Let {uε } be a bounded sequence in L2 (Ω ). Then the following propositions are equivalent: (i) Tε (uε ) weakly converges to u0 in L2 (Ω × Y ). (ii) uε two-scale converges to u0 . Similarly to the two-scale convergence technique, the following properties hold: Proposition 4. (i) Let uε be a bounded sequence in L2 (Ω ). Then there exists u0 ∈ L2 (Ω × Y ) such that, up to a subsequence, Tε (uε ) weakly converges to u0 in L2 (Ω × Y ). (ii) Let uε be a bounded sequence in H 1 (Ω ), which weakly converges to a limit u0 ∈ H 1 (Ω ). Then Tε (uε ) weakly converges to u0 1 in L2 (Ω × Y ) and there exists a function u1 ∈ L2 (Ω ; Hper (Y )/R) such that, up to a subsequence, Tε (∇ uε ) weakly converges to ∇x u0 + ∇y u1 in L2 (Ω × Y ).

To conclude this brief introduction, we have the following (straightforward) property: Proposition 5. Let u ∈ L2 (Ω ). If a sequence uε ∈ L2 (Ω ) strongly converges to u in L2 (Ω ), then Tε (uε ) strongly converges to u in L2 (Ω × Y ).

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

2137

Fig. 2. Rectangular dam with layers.

2.2. Estimates and convergence results Using the periodic unfolding method, we establish the behavior of the fluid flow when ε goes to 0. From now on let us consider a dam with layers (see Fig. 2). Let us suppose that: Assumption 2. The permeability kε of the porous medium can be written as: kε (x) :=

K11−α



γ

x,

X1 (x)

ε



γ



α

x,

K2

X2 (x)

 (3)

ε

with α ∈ [0, 1]. Here, the mapping Fγ : x → Xγ (x) denotes the rotation of angle γ , that is:

 Fγ :

γ

X1 (x) = cos γ x1 − sin γ x2 , γ

X2 (x) = sin γ x1 + cos γ x2 .

We should notice that γ = 0 (resp. γ = π /2) may lead to the consideration of horizontal and vertical layers. In the forthcoming section of numerical results, despite examples of these two cases, also an example of oblique layers with γ = π/3 are presented. Definition 3. The homogenized dam problem is defined by:

 Find (p0 , Θ1 , Θ2 ) ∈ V0∩ Wα × L∞ (Ω ) × L∞ (Ω ) such that   (A0 · ∇ p0 ) · ∇φ − (B 0 (Θ1 , Θ2 ) · g) · ∇φ ≤ 0, ∀φ ∈ V+ ∩ W0 , (A⋆ )  Ω  Ω p0 ≥ 0, Θi ∈ H (p0 ), (i = 1, 2), a.e. in Ω , with the homogenized permeability matrix

A0 =



κ1⋆

0



κ2⋆

0

− (κ1⋆ − κ2⋆ ) sin γ

sin γ cos γ



cos γ − sin γ



,

and the homogenized gravity term:

B (Θ1 , Θ2 ) = 0



Θ1 κ1⋆ 0



0

Θ2 κ2⋆

   sin γ ⋆ ⋆ − Θ1 κ1 − Θ2 κ2 sin γ cos γ

cos γ − sin γ



.

The homogenized coefficients are defined by ⋆

κ1 =

K2α

Y

K1−1+α

Y

,

⋆

κ2 =

K11−α K2−α

Y

Y

.

Theorem 1. There exists (p0 , θ0 ) ∈ V0 ∩ Wα × L∞ (Ω × Y ) such that

• p0 ≥ 0 and θ0 ∈ H (p0 ) a.e.; • (p0 , θ0 ) is the two-scale limit of (pε , θε ), up to a subsequence of ε . Note that the subsequence is still denoted by ε for the sake of simplicity.

2138

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

Defining the following saturations

 Θ1 (x) =

Y θ0 (x, ·)K2α (x, ·) Y

K2α (x, ·)

,

Θ2 (x) =

 Y θ0 (x, ·)K11−α (x, ·) Y

K11−α (x, ·)

,

(4)

then (p0 , Θ1 , Θ2 ) is a solution of problem (A⋆ ). Proof. The convergence result is proved using different steps: 1. 2. 3. 4. 5. 6. 7.

Formulation of the dam problem in the layer-structured coordinates. Estimates, convergence results and properties of the limits. Macroscopic and microscopic decompositions. Definition of the cell problems. Definition of the homogenized problem. Computation of the homogenized coefficients. Formulation of the homogenized problem in the initial coordinates.

First let us establish that the solution is bounded with respect to ε . 1. Formulation of the dam problem in the layer-structured coordinates. Using the change of coordinates x → Fγ (x), we briefly introduce the following notation

 f ( x) = f ◦ Fγ−1 ( x). Then the dam problem may be rewritten as

 ∞    Find ( pε ,  θε ) ∈  V0 ∩   Wα × L (Ω ) such that:   0 , kε ∇ pε · ∇ φ−  θε kε gγ · ∇ φ ≤ 0, ∀ φ ∈ V+ ∩ W ( Aε )    Ω Ω   , pε ≥ 0,  θε ∈ H ( pε ), a.e. in Ω with the following functional spaces:

    ), v|Γ = 0 , V0 = v ∈ H 1 (Ω 0   α = v ∈ H 1 (Ω  ), v|Γ =  W α , a

    ), v|Γ ≥ 0 , V+ = v ∈ H 1 ( Ω 0   0 = v ∈ H 1 (Ω  ), v|Γ = 0 . W a

Let us point out that the change of coordinates only has affected the effective gravity vector, as g has been changed into γ



g :=

cos γ sin γ

− sin γ cos γ



· g.

In the sequel, for the sake of simplicity, overscripts will be dropped. 2. Estimates, convergence results and properties of the limits. Since 0 ≤ θε ≤ 1, θε is bounded in L∞ (Ω ) and in L2 (Ω ). Moreover, using pε −  α (with  α a regular function such that pε −  α ∈ V+ ∩ W0 ) as a test function, Poincaré–Friedrichs and Cauchy–Schwarz inequalities, we prove that pε is bounded in H 1 (Ω ). Thus, as a consequence of the previous estimates (see Propositions 3 and 4, or Proposition 1.14 in [29], Theorem 1 13 in [30]), there exists (p0 , p1 ) ∈ V0 ∩ Wα × L2 (Ω ; Hper (Y )/R) and θ0 ∈ L2 (Ω × Y ) such that, up to a subsequence, the

following convergences hold in L2 (Ω × Y ):

Tε (pε ) → p0 ,

Tε (∇ pε ) ⇀ ∇x p0 + ∇y p1 ,

Tε (θε ) ⇀ θ0 .

Let us prove that p0 ≥ 0,

θ0 ∈ H (p0 ), a.e.

(5)

As pε ≥ 0 and 0 ≤ θε ≤ 1 a.e., using the definition of the unfolding operator, on has: Tε (pε ) ≥ 0 and 0 ≤ Tε (θε ) ≤ 1 a.e. Passing to the limit on ε yields p0 ≥ 0,

0 ≤ θ0 ≤ 1,

a.e. in Ω × Y .

Applying the unfolding operator to the equality pε (1 −θε ) = 0 and passing to the limit, we get: p0 (1 − θ0 ) = 0 in L1 (Ω × Y ). Since p0 ≥ 0 and (1 − θ0 ) ≥ 0 a.e., and recalling that pε (resp. θε ) strongly (resp. weakly) converges to p0 (resp. θ0 ): p0 (1 − θ0 ) = 0,

a.e. in Ω × Y .

Thus, we have either (p0 > 0 and θ0 = 1) or (p0 = 0 and 0 ≤ θ0 ≤ 1), so that θ0 ∈ H (p0 ) a.e., which show the thesis.

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

2139

3. Macroscopic and microscopic decompositions. The limit functions p0 , p1 and θ0 satisfy the following equations:

• Macroscopic equation. For every φ in V+ ∩ W0 ,         K (x, y) ∇x p0 (x) + ∇y p1 (x, y) dy · ∇φ(x) dx − θ0 (x, y)K (x, y) dy gγ · ∇φ(x) dx ≤ 0. Ω

Ω

Y

Y

1 • Microscopic equation. For a.e. x ∈ Ω , for every ψ ∈ Hper (Y ),     K (x, y) ∇x p0 (x) + ∇y p1 (x, y) · ∇ψ(y) dy − θ0 (x, y)K (x, y) gγ · ∇ψ(y) dy = 0. Y

Y

Considering problem (Aε ), we have, using Proposition 2:

 Ω ×Y

Tε (kε ) Tε (∇ pε ) · Tε (∇ϕε ) −

 Ω ×Y

Tε (θε )Tε (kε ) gγ · Tε (∇ϕε ) ≤ 0,

for every ϕε ∈ V+ ∩ W0 . Taking a test function x → φ(x), with φ ∈ V+ ∩ W0 and using the convergence results established in the previous step gives the macroscopic equation by passing to the limit with respect to ε . Now taking the test functions 1 x → ± εφ(x) ψ(x/ε), with φ ∈ Cc∞ (Ω ) and ψ ∈ Hper (Y ), gives the microscopic one. 4. Definition of the cell problems. Let us define the local problems, respectively denoted (Ni⋆ ) and (Ni0 ):

 ⋆ 2 1 Find  Wi (i = 1, 2) in L (Ω ; Hper(Y )/R) such that, for a.e. x ∈ Ω , ⋆ ∂ψ (Ni ) 1 , ∀ψ ∈ Hper (Y ).  K (x, ·) ∇y Wi⋆ (x, ·) · ∇y ψ = K (x, ·) ∂ y i Y Y  1  Find Wi0 (i = 1, 2) in L2 (Ω ; Hper (Y )/R) such that, for a.e. x ∈ Ω , 0 ∂ψ (Ni ) 1 0 , ∀ψ ∈ Hper (Y ).  K (x, ·) ∇y Wi (x, ·) · ∇y ψ = θ0 (x, ·)K (x, ·) ∂ yi Y Y As a consequence of Lax–Milgram theorem, local problem (Ni⋆ ) (resp. (Ni0 )) admits a unique solution Wi⋆ (resp. Wi0 ) in 1 L (Ω ; Hper (Y )/R). 5. Definition of the homogenized problem. 1 From the local problems, we easily obtain in L2 (Ω ; Hper (Y )/R): 2

W1⋆ (x, y) p1 (x, y) = − W2⋆ (x, y)





· ∇x p0 (x) +



W10 (x, y) W20 (x, y)



· gγ .

The homogenized formulation follows by substituting the previous expression of p1 in the macroscopic equation:

 Ω

A0 · ∇ p0 ∇φ −



B 0 · gγ · ∇φ ≤ 0,

 Ω



∀φ ∈ V+ ∩ W0 ,

with

A0ij

:= K δij − K

∂ Wj⋆ ∂ yi

Y

,

Bij0

:= θ0 K δij − K

∂ Wj0 ∂ yi

Y

.

6. Computation of the homogenized coefficients. Let us recall the layer-structured assumption: K (x, y) := K11−α (x, y1 ) K2α (x, y2 ). Let us compute the permeability matrix.

• First, considering problem (N1⋆ ), one gets:   ∂ψ 1 K ∇y W1⋆ · ∇y ψ = K , ∀ψ ∈ Hper (Y ). ∂ y1 Y Y Next, using a test function only depending on y1 , one has: Y2

 K Y1

∂ W1⋆ ψ′ = ∂ y1



Y2

K ψ ′, Y1

1 ∀ψ ∈ Hper (Y1 ).

2140

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

Then, for a.e. x ∈ Ω K (x, ·)

∂ W1⋆ (x, ·) ∂ y1

Y2

Y2

= K (x, ·) − C1⋆ (x),

(6)

where C1⋆ (x) is an additive constant with respect to y. Using the layer-structured assumption, dividing by K11−α , averaging with respect to the y2 variable and using the Y periodicity of W1⋆ , leads to the following equality: Y

Y

K2α (x, ·) − C1⋆ (x) K1−1+α (x, ·) = 0.

(7)

Next, from the definition of A0 (see the previous step) and Eq. (6), we deduce that C1⋆ (x) = A011 (x) so that we get from Eq. (7):

A011 (x) =

K2α (x, ·)

Y Y

K1−1+α (x, ·)

.

A straightforward adaptation of the previous computation allows us to conclude that Y

A022

(x) =

K11−α (x, ·) K2−α (x, ·)

Y

.

• Second, considering problem (N1⋆ ), one gets:   ∂ψ 1 ⋆ , ∀ψ ∈ Hper (Y ). K ∇y W1 · ∇y ψ = K ∂ y1 Y Y Next, using a test function only depending on y2 , one has: Y1

 K Y2

∂ W1⋆ ψ ′ = 0, ∂ y2

1 ∀ψ ∈ Hper (Y2 ).

Then, for a.e. x ∈ Ω K (x, ·)

∂ W1⋆ (x, ·) ∂ y2

Y1

= −C2⋆ (x),

(8)

where C2⋆ (x) is an additive constant with respect to y. Using the layer-structured assumption, dividing by K2α , averaging with respect to the y1 variable and using the Y periodicity of W1⋆ , leads to the following equality: C2⋆ (x) = 0.

(9)

Next, from the definition of A0 (see the previous step) and Eq. (8), we deduce that C2⋆ (x) = A021 (x) so that we get from Eq. (9):

A021 (x) = 0. In the same way, we have

A012 (x) = 0. Let us compute the right-hand side matrix.

• First, considering problem (N10 ), one gets:   ∂ψ 1 , ∀ψ ∈ Hper (Y ). K ∇y W10 · ∇y ψ = θ0 K ∂ y1 Y Y Next, using a test function only depending on y1 , one has: Y2

 K Y1

∂ W10 ψ′ = ∂ y1



Y

θ0 K 2 ψ ′ ,

1 ∀ψ ∈ Hper (Y1 ).

Y1

Then, for a.e. x ∈ Ω Y2

K (x, ·)

∂ W10 (x, ·) ∂ y1

Y2

= θ0 (x, ·)K (x, ·) − C10 (x),

(10)

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

2141

where C10 (x) is an additive constant with respect to y. Using the layer-structured assumption, dividing by K11−α , averaging with respect to the y2 variable and using the Y periodicity of W10 , leads to the following equality: Y

Y

θ0 (x, ·)K2α (x, ·) − C10 (x) K1−1+α (x, ·) = 0.

(11)

0 Next, from the definition of B 0 (see the previous step) and Eq. (10), we deduce that C10 (x) = B11 (x) so that we get from Eq. (11): Y

0 B11 (x) =

θ0 (x, ·)K2α (x, ·)

.

Y

K1−1+α (x, ·)

A simple renormalization allows us to define Y

0 B11

(x) = Θ1 (x)

Y

K2α (x, ·)

K1−1+α (x, ·)

Y

,

with Θ1 (x) :=

θ0 (x, ·)K2α (x, ·) K2α (x, ·)

Y

.

A straightforward adaptation of the previous computation allows us to conclude that 0 B22

(x) = Θ2 (x)

K11−α (x, ·)

Y

Y

K2−α (x, ·)

,

with Θ2 :=

θ0 K11−α K11−α

Y

Y

.

• Second, considering problem (N10 ), one gets:   ∂ψ 1 K ∇y W10 · ∇y ψ = θ0 K , ∀ψ ∈ Hper (Y ). ∂ y 1 Y Y Next, using a test function only depending on y2 , one has:



Y1

∂ W10 ψ ′ = 0, K ∂ y2 Y2

1 ∀ψ ∈ Hper (Y2 ).

Then, for a.e. x ∈ Ω Y1

∂ W10 K (x, ·) (x, ·) ∂ y2

= −C20 (x),

(12)

where C20 (x) is an additive constant with respect to y. Using the layer-structured assumption, dividing by K2α , averaging with respect to the y1 variable and using the Y periodicity of W10 , leads to the following equality: C20 (x) = 0.

(13) 0

Next, from the definition of B (see the previous step) and Eq. (12), we deduce that Eq. (13):

C20

(x) =

0 B21

(x) so that we get from

0 B21 (x) = 0.

In the same way, we have 0 B12 (x) = 0.

Finally, as 0 ≤ θ0 ≤ 1 and K2 > 0, one gets Y

0 ≤ θ0 (x, ·)K2α (x, ·) ≤ K2α (x, ·)

Y

so that 0 ≤ Θ1 ≤ 1. Since θ0 ∈ H (p0 ), Θ1 = 1 if θ0 = 1 and Θ1 < 1 if θ0 < 1, we finally obtain

Θ1 ∈ H (p0 ). In the same way,

Θ2 ∈ H (p0 ). 7. Formulation of the homogenized problem in the initial coordinates. Synthesizing the previous steps (and reintroducing the notations related to the layer-structured coordinates), the homogenized formulation is defined by:

 α × L∞ (Ω 1 , Θ 2 ) ∈   ) × L∞ ( Ω  ) such that Find ( p0 , Θ V0∩ W   0 0 γ 0 , 1 , Θ 2 ) · g ) · ∇  · ∇  (Θ (A p0 ) · ∇ φ − (B φ ≤ 0, ∀ φ ∈ V+ ∩ W    Ω  Ω i ∈ H ( ,  p0 ≥ 0, Θ p0 ), (i = 1, 2), a.e. in Ω

2142

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

where the homogenized permeability and gravity matrices take the following form

 ⋆  κ1 0  A = 0

0



 κ2⋆

 1  Θ κ1⋆ 0    B (Θ1 , Θ2 ) =

,

0



, 2  Θ κ2⋆

0

and with the following homogenized coefficients

 κ1⋆ ( x) :=

 K α ( x, ·)

Y

Y

, Y  K1−1+α ( x, ·) 2

 κ2⋆ ( x) :=

 K 1−α ( x, ·)

.

1

 K −α ( x, ·)

Y

2

Besides, the link between the homogenized saturation functions and the microscopic (two-scale limit) saturation function is given by

1 ( Θ x) =

 Y  x, ·) K2α ( x, ·) θ0 ( Y

 K α ( x, ·)

,

2 ( Θ x) =

 Y  x, ·) K11−α ( x, ·) θ0 (

2

Y

 K 1−α ( x, ·)

.

1

Using the change of coordinates  x →

Fγ−1 ( x)

and, consequently, the following notation

f (x) =  f ◦ Fγ (x), we are now able to rewrite the homogenized problem in the initial coordinates, which leads us to conclude the proof after simple computations.  Theorem 2. Problem (A⋆ ) admits a so-called ‘‘isotropic solution’’, i.e. (A⋆ ) admits a solution (p, Θ1 , Θ2 ) such that Θ1 = Θ2 . Proof. Existence of an isotropic solution relies on the homogenization of the penalized problem. Let us recall the definition of the penalized problem (with oscillating permeability):

 Find pηε ∈ V0 ∩ Wα such that:   k ∇ pηε · ∇φ − Hη (pηε )k g · ∇φ ≤ 0, (Aεη )  Ω Ω  η pε ≥ 0, a.e. in Ω .

∀φ ∈ V+ ∩ W0 , η

η

Following the same guideline as in the initial problem, we prove that there exists p0 ∈ V0 ∩ Wα , p0 ≥ 0 a.e., which is the η two-scale limit of pηε up to a subsequence. Moreover, the limit p0 satisfies the following homogenized penalized problem:

 η Find p0 ∈ V0 ∩ Wα such   that   η η η 0 ⋆ (A · ∇ p0 ) · ∇φ − (B 0 (Hη (p0 ), Hη (p0 )) · g) · ∇φ ≤ 0, (Aη )  Ω Ω   η p0 ≥ 0, a.e. in Ω .

∀φ ∈ V+ ∩ W0 ,

η

Then, using a suitable test-function, we prove that the sequence {p0 }η is bounded, in the H 1 -norm, by a constant which η does not depend on η. Besides, due to L∞ -boundedness of Hη (p0 ), we obtain the following convergence results (up to a ∞ subsequence): there exists p0 ∈ V0 ∩ Wα , Θ ∈ L (Ω ) such that η

• p0 ⇀ p0 in H 1 (Ω ), η • Hη (p0 ) ⇀ Θ in L∞ (Ω ) − ⋆. Passing to the limit with respect to. η in (A⋆η ), we obtain

 Ω

(A0 · ∇ p0 ) · ∇φ −

 Ω

(B 0 (Θ , Θ ) · g) · ∇φ ≤ 0,

∀φ ∈ V+ ∩ W0 . η

Besides, nonnegativity of p0 is inherited from the properties of p0 . So are the bounds of Θ (i.e. Θ ∈ [0, 1]) which follow from η η the properties of the penalized Heaviside function. Finally, one can prove that p0 (1 − Θ ) = 0 by checking that p0 (1 − Hη (p0 )) tends to 0. This allows us to conclude that Θ ∈ H (p0 ).  Corollary 3 (Dam Problem with Horizontal and Vertical Layers). The homogenized dam problem, for g := (0, −g ) and γ := 0, is defined by:

 ∞ Find (p, Θ ) ∈ V0 ∩ W α × L (Ω ) such that   (A0 · ∇ p) · ∇φ − (Θ κ2⋆ ) g · ∇φ ≤ 0, (A⋆ )  Ω  Ω p ≥ 0, Θ ∈ H (p), a.e.

∀φ ∈ V+ ∩ W0 ,

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

2143

Moreover problem (A⋆ ) admits (p0 , Θ ) as a solution, with

 Θ=

θ0 K11−α K11−α

Y

Y

,

where (p0 , θ0 ) is the two-scale limit of (pε , θε ) (solution of problem (Aε )). Remark 2. For g := (0, −g ) and γ := 0, uniqueness of the solution (p, Θ ) can be obtained under additional assumptions: with the terminology used by Chipot [31], ‘‘if the shape of Ω is such that no pool appears’’ and if the homogenized coefficients are constant (see [32]), then the solution is unique (it is quite easy to make a change of coordinates so that we get an homogeneous dam problem) and the homogenized saturation Θ is a characteristic function. Remark 3. For g := (0, −g ), γ := 0 and α := 1, then the solution of the homogenized problem is the weak limit (in L2 ) of the solution of the initial problem (up to a subsequence). We emphasize that the homogenization of the dam problem with leaky boundary conditions can be readily adapted from the previous guideline. For this purpose, let us consider the initial problem with oscillating coefficients:

 ∞ Find (pε , θε ) ∈ V0 ×   L (Ω ) such that:   kε ∇ pε · ∇φ − θε kε g · ∇φ ≤ β(·, α − pε ) φ, (Bε )  Ω Γa  Ω pε ≥ 0, θε ∈ H (pε ), a.e. in Ω .

∀φ ∈ V+ ,

Definition 4. The homogenized dam problem with leaky boundary conditions is defined by:

 Find (p0 , Θ1 , Θ2 ) ∈ V0× L∞ (Ω ) × L∞ (Ω ) such that    0 0 ⋆ ( A · ∇ p ) · ∇φ − ( B ( Θ , Θ ) · g ) · ∇φ ≤ β(·, α − p0 ), (B ) 0 1 2  Ω Γa  Ω p0 ≥ 0, Θi ∈ H (p0 ), (i = 1, 2), a.e. in Ω ,

∀φ ∈ V+ ,

where the homogenized coefficients are given by Definition 3. Theorem 4. There exists (p0 , θ0 ) ∈ V0 × L∞ (Ω × Y ) such that

• p0 ≥ 0 and θ0 ∈ H (p0 ) a.e.; • (p0 , θ0 ) is the two-scale limit of (pε , θε ), up to a subsequence of ε . Note that the subsequence is still denoted by ε for the sake of simplicity. Defining the saturation functions (Θi )i=1,2 by Eq. (4), then (p0 , Θ1 , Θ2 ) is a solution of problem (B⋆ ). Besides, (B⋆ ) admits an isotropic solution. Proof. The guideline of the proof is a straightforward adaptation of the previous results dealing with Dirichlet boundary conditions. When leaky boundary conditions are taken into account, we show that pε is bounded in H 1 (Ω ), by adapting the computations developed in [22] for the regularized problem. The, using classical compactness arguments [33], there exists p ∈ V0 ⊂ H 1 (Ω ) such that, up to a subsequence, pε ⇀ p in H 1 (Ω ), pε → p in L2 (Ω ) and a.e. on Ω , and γa (pε ) → γa (p) in L2 (Γa ), where γa denotes the trace operator on Γa . Then, the previous analysis conducted with Dirichlet boundary conditions can be adapted and only the behavior of the boundary term needs to be further detailed. More precisely, due to uniform Lipschitz continuity of β with respect to the second variable, Cauchy–Schwartz inequality and strong convergence of γa (pε ) in L2 (Γa ), we get: for any φ ∈ V+ ,

 lim

ε→0 Γ a

β(·, α − pε )φ =

which concludes the proof.

 Γa

β(·, α − p0 )φ,



3. Numerical methods In this section we describe the numerical methods to solve the different parameter dependent problems and the corresponding homogenized ones. These methods are mainly an adaption to the here considered variable and nondiagonal permeability matrix setting of the numerical techniques proposed in Bermúdez and Durany [5] for the case of constant permeability, which are based on the application of characteristics methods to steady state convection–diffusion equations with a nonlinear convection term, the nonlinearity being motivated by the presence of the Heaviside multivalued operator. These kind of techniques have already been used in the related free boundary problem which is defined by the Elrod–Adams model for cavitation in lubrication (see [34–36], the last work corresponding to an homogenization study related to lubrication of surfaces with periodic roughness).

2144

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

3.1. Characteristics method for the convective term In order to describe the application of the characteristics method for the particular convective term in the convection–diffusion equation, we first take into account that div g = 0 and consider the following equivalent form of Eq. (1): g · ∇ (θ k) − div (k∇ p) = 0,

p ≥ 0, θ ∈ H (p).

(14)

Next, in order to apply the method of characteristics, we consider the artificial time dependent velocity field u(t , x) = g and an artificial time dependence in all involved functions, by introducing the notation φ(t , x) = φ(x). Then, we can introduce the following material (or total derivative) operator associated to the velocity field u: D Dt

=

∂ + u · ∇. ∂t

(15)

So, we can write Eq. (16) in the form D  Dt

   θ k − div k∇ p = 0,

p ≥ 0, θ ∈ H (p).

(16)

Next, by approximating the total derivative by the characteristics method, the previous equation can be approximated by

  θ n+1 kn+1 − (θ n kn ) ◦ κ n − div kn+1 ∇ pn+1 = 0, 1t

pn+1 ≥ 0, θ n+1 ∈ H (pn+1 ),

where the index n + 1 denotes the approximation at the artificial time t n+1 of the introduced time dependent functions and 1t denotes an artificial time step. Moreover, the function κ n is defined at each spatial point of the domain by κ n (x) = κ(t n+1 , x; t n ), that denotes the position at time t n of the point placed in x at time t n+1 and moving along the integral path (characteristic curve) defined by the velocity field u, so that κ n (x) can be obtained from the solution of the final value ODE problem: dκ dτ

  (t n+1 , x; τ ) = u τ , κ(t n+1 , x; τ ) ,

κ(t n+1 , x; t n+1 ) = x.

For the particular expression of the velocity field u, for x = (x1 , x2 ) we can easily compute

κ n (x1 , x2 ) = (x1 , x2 + g 1t ).

(17)

Moreover, removing the artificial time dependence of the permeability k, we get

  θ n+1 k − (θ n k) ◦ κ n − div k∇ pn+1 = 0, 1t

pn+1 ≥ 0, θ n+1 ∈ H (pn+1 ).

(18)

In order to pose an approximated problem (which is actually dependent on the parameter 1t), Eq. (18) is completed with the boundary conditions in (2) now written in terms of pn+1 and θ n+1 , instead of p and θ . Next, the variational formulation for the approximated problem defined by (18) can be written as follows: Find (pn+1 , θ n+1 ) ∈ V− × L∞ (Ω ) such that pn|Γ+a1 = α, θ n+1 ∈ H (pn+1 ) and

1t

 Ω

k ∇ pn+1 · ∇(φ − pn+1 ) +

+ 1t

 Γ0 ∪Γ



θ n+1 k (φ − pn+1 )   n  n+1 n +1 θ k n · e2 (φ − p ) ≥ (θ k) ◦ κ n (φ − pn+1 ), Ω

Ω

for all φ ∈ V− such that φ|Γa = 0, where V− = {v ∈ H 1 (Ω ), v|Γ ≤ 0}. Next, by using the indicatrix of the convex set V− , 0 denoted by IV− , we get the more convenient formulation: Find (pn+1 , θ n+1 ) ∈ V− × L∞ (Ω ) such that pn|Γ+a1 = α, θ n+1 ∈ H (pn+1 ) and

1t

 k ∇p Ω

 ≥ Ω

n+1

· ∇(φ − p

n+1

)+

 Ω

θ

n+1

k (φ − p

n+1

) + 1t

 Γ0 ∪Γ

θ n+1 k n · e2 (φ − pn+1 ) + IV− (φ) − IV− (pn+1 )

 n  (θ k) ◦ κ n (φ − pn+1 ),

for all φ ∈ H 1 (Ω ) such that φ|Γa = 0. Next, by using the definition of subdifferential operator [37], we get the following formulation:

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

2145

Find (pn+1 , θ n+1 ) ∈ V0 ∩ Wα × L∞ (Ω ) such that:

1t



k∇ pn+1 · ∇φ +



Ω

Ω

θ n+1 k φ + 1t

 Γ0 ∪Γ

θ n+1 k n · e2 φ + 1t

 Γ0

qn+1 φ =

 Ω

 n  (θ k) ◦ κ n φ,

(19)

for all φ ∈ H 1 (Ω ) such that φ|Γa = 0 and qn+1 ∈ ∂ IV− (pn+1 ),

(20)

θ

(21)

n+1

∈ H (p

n +1

),

where ∂ IV− (pn+1 ) denotes the subdifferential operator of IV− at the point pn+1 . 3.2. Duality method for the nonlinear terms In this section we describe the application of a duality algorithm proposed in [38] to treat numerically the nonlinearities appearing in Eqs. (19)–(21). For this purpose, we first notice that the nonlinear relations defined by Eqs. (20) and (21) are expressed in terms of the multivalued maximal monotone operators ∂ IV− and H, respectively. We address the reader to [37] for further details about the definitions and results concerning maximal monotone operators we use in this paper. Let G be a maximal monotone operator and let ω and λ be two nonnegative parameters satisfying the condition λω < 1, then it can be proved that the resolvent operator of G, Jλω = ((1 − λω)I + λG)−1 , is well defined. Furthermore, the corresponding Yosida approximation of operator G − ωI is defined by Gω λ =

I − Jλω

λ

and it is a Lipschitz operator with constant λ−1 . In this paper, we propose an algorithm based on the following result (see [38], for the details): Lemma 6. Let G be a maximal monotone operator in a Hilbert space V . Then the following conditions are equivalent: 1. u ∈ G(y) − ωy. 2. u = Gω λ (y + λu). Next, in terms of the nonnegative parameters ω1 and ω2 , we introduce the following new unknowns (multipliers):

α n+1 ∈ ∂ IV− (pn+1 ) − ω1 pn+1 ,

(22)

β n+1 ∈ H (pn+1 ) − ω2 pn+1 .

(23)

So, if we apply Lemma 6 to the maximal monotone operators ∂ IV− and H, then problem (19)–(21) can be equivalently written in the form:



k∇ pn+1 · ∇φ + Ω

=

1

1t



 Ω

ω2 1t

(θ n k) ◦ κ



 n

Ω

kpn+1 φ +

φ−

1



1t

Ω

 Γ0 ∪Γ

θ n+1 k n · e2 φ + ω1

β n +1 k φ −

 Γ0

α n+1 φ,

 Γ0

pn+1 φ (24)

for all φ ∈ H 1 (Ω ) such that φ|Γa = 0 and ω

α n+1 = (∂ IV− )λ 1 (pn+1 + λ1 α n+1 ), β

n+1

=

ω Hλ 2 (pn+1

+ λ2 β

n+1

).

(25) (26)

Notice that by using Lemma 6 we have replaced a nonlinear multivalued relation by a nonlinear univalued one. The algorithm allows us to solve the remaining nonlinearities in problem (24)–(26) by means of a fixed-point iteration that alternatively solves (24) for given values of the multipliers α n+1 and β n+1 , and updates these multipliers from the last computed value of pn+1 by using (25) and (26). The fixed-point iteration is initialized by prescribing the initial values of the multiplier at iteration n + 1 to the ones obtained after convergence at the previous iteration n. For convergence purposes, the values of λi and ωi are chosen so that λi ωi = 0.5, i = 1, 2, this is the reason why in next section we just indicate the value of parameters ω ω ωi . Also, the functions (∂ IV− )λ 1 and Hλ 2 can be explicitly computed. For the spatial discretization of the weak formulation we consider a triangular mesh τh of the domain Ω , so that the space H 1 (Ω ) is approximated by the following piecewise linear Lagrange finite elements space:

¯ ) : φh |T ∈ P1 , ∀T ∈ τh }, Vh = {φh ∈ C 0 (Ω

(27)

¯ ) is the space of continuous functions on Ω ¯ and P1 denotes the space of polynomials of degree less or equal where C 0 (Ω than one.

2146

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151 Table 1 Parameters of the model and the numerical methods that have been used in the simulation results. Figs.

3–5

6–8

9–11

Domain Width L Height H

Rectangular 1.00 2.00

id. id. id.

id. id. id.

Gravity g

(0, −1)

id.

id.

1.90 0.10

id. id.

id. id.

251 601 150 851 300 000

601 251 150 851 300 000

381 381 145 161 288 800

2.0 · 10−3 1.0 · 10−6 0.5 0.5 1.0 · 10−5 1.0 · 10−6

id. id. id. id. id. id.

id. id. id. id. id. id.

Boundary conditions Water height (left wall) Water height (right wall) Mesh Number of vertices on horizontal walls Number of vertices on lateral walls Number of vertices Number of triangles Algorithm parameter Time step 1t (characteristics) Tolerance in saturation (characteristics) Parameter ω1 (duality) Parameter ω2 (duality) Tolerance in multipliers (duality) Tolerance in pressure (duality)

3.3. Algorithm implementation For simplicity, in previous sections concerning the description of the numerical techniques we have explained their application to the standard dam problem defined by Eq. (1). In practice, the set of numerical techniques has been applied also to the ε -dependent problems associated to the stratified porous medium permeability defined by Eq. (3) and the corresponding homogenized dam problem (A⋆ ). Taking into account the schemes of the previously described methods, the global algorithm for solving any of the problems (the parameter dependent and the homogenized ones) contains two loops: the inner one corresponding to the duality method for the nonlinearities, and the outer one corresponding to the use of characteristics method to deal with the convection term. The outer loop can also be understood as intended to converge from the solution of the artificial transient problem to the initial posed steady-state one. Concerning the characteristics loop we notice that the artificial time step has been chosen such that g 1t = 1y, in order to guarantee that κ n (vi ) also belongs to the mesh for each mesh node vi , thus avoiding the need of additional interpolations to compute ((k θ n ) ◦ κ n )(vi ). Furthermore, as expression (17) shows that κ n (vi ) actually does not depend on n, these computations can be performed just once outside the characteristics loop. The same occurs with the matrix computation and the factorization of the linear system to be solved at each step of the inner duality method loop. Therefore, just the second members of the linear systems have to be updated in the multiplier and characteristics loops accordingly. Concerning the duality method loop, the values of α n+1 and β n+1 in the mesh nodes are updated and provided to the linear problem at each fixed point iteration, starting from their initial values taken from the last ones of the previous characteristics loop iteration. 4. Numerical results In this section, the numerical simulation of different layered-structured rectangular dams is performed to illustrate the theoretical results of convergence that have been stated in Section 2, as well as the good performance of the described numerical techniques. For this purpose, the parameters of the model and numerical methods appearing in Table 1 have been considered. Besides the numerical results here presented, appropriate numerical tests have been performed to validate the convergence as soon as 1t decreases and the spatial mesh becomes finer. 4.1. Simulation results: direct problems and homogenized problems In particular we compute the solution of the dam problem with horizontal, vertical and oblique layers and compare the corresponding solution to the homogenized one. All the simulations have been performed with data provided by Table 1. In next paragraphs the obtained results are described and discussed:

• Horizontal layers. This test example corresponds to the particular case x  2 , kε (x) := K2 ε

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

2147

1

1.8

0.9

1.6

0.8 1.4 0.7 1.2

0.6

1

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

i)

ii) Fig. 3. (i) Pressure and (ii) saturation for horizontal layers and ε

0 −1

= 10. 1

1.8

0.9

1.6

0.8 1.4

0.7

1.2

0.6

1

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

i)

ii) Fig. 4. (i) Pressure and (ii) saturation for horizontal layers and ε

0 −1

= 30. 1

1.8

0.9

1.6

0.8 1.4 0.7 1.2

0.6

1

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

i)

ii)

0

Fig. 5. (i) Pressure and (ii) saturation for horizontal layers in the homogenized case.

i.e. γ = 0 and α = 1, with different values of ε . More precisely, we have defined K2 (y) =



10, 1,

if 0.0 ≤ y < 0.5, if 0.5 ≤ y < 1.0.

Computing the homogenized solution leads us to consider the following permeability matrix and gravity term, with g = g (0, −1): 0

A ≃



5.50 0.00

0.00 , 1.82



⋆

Θ κ2 g ≃ −Θ

0.00 . 1.82





Figs. 3–5 show the distribution of the pressure (isobaric lines) and saturation in the dam domain for different values of the period, namely ε = 1/10, 1/30 and in the homogenized case.

2148

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

1

1.8

0.9

1.6

0.8 1.4 0.7 1.2

0.6

1

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

i)

ii)

0

Fig. 6. (i) Pressure and (ii) saturation for vertical layers and ε −1 = 10.

1

1.8

0.9

1.6

0.8 1.4

0.7

1.2

0.6

1

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

i)

ii) Fig. 7. (i) Pressure and (ii) saturation for vertical layers and ε

0 −1

= 30.

1

1.8

0.9

1.6

0.8 1.4

0.7

1.2

0.6

1

i)

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

ii) Fig. 8. (i) Pressure and (ii) saturation for vertical layers in the homogenized case.

• Vertical layers. This test example corresponds to the particular case x  1 kε (x) := K1 , ε with, for instance, γ = 0 and α = 0. More precisely, we have defined K1 (y) =



10, 1,

if 0.0 ≤ y < 0.5, if 0.5 ≤ y < 1.0.

0

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

2149

1

1.8

0.9

1.6

0.8 1.4 0.7 1.2

0.6

1

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

i)

ii)

0

Fig. 9. (i) Pressure and (ii) saturation for oblique layers and ε −1 = 10.

1

1.8

0.9 1.6 0.8 1.4

0.7

1.2

0.6

1

0.5

0.8

0.4

0.6

0.3

0.4

0.2 0.1

0.2

i)

ii)

0

Fig. 10. (i) Pressure and (ii) saturation for oblique layers and ε −1 = 30.

1

1.8

0.9

1.6

0.8 1.4

0.7

1.2

0.6

1

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

i)

ii)

0

Fig. 11. (i) Pressure and (ii) saturation for oblique layers in the homogenized case.

Computing the homogenized solution leads us to consider the following permeability matrix and gravity term, with g = g (0, −1): 0

A ≃



1.82 0.00

0.00 , 5.50



⋆

Θ κ2 g ≃ −Θ



0.00 . 5.50



Figs. 6–8 show the distribution of the pressure (isobaric lines) and saturation in the dam domain for different values of the period, namely ε = 1/10, 1/30 and in the homogenized case. • Oblique layers. This test example corresponds to the particular case: kε (x) := K2



sin(−π /3)x1 + cos(−π /3)x2

ε



.

2150

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

Fig. 12. Order of convergence of the pressure in the discrete L2 norm for (A) horizontal layers, (B) vertical layers and (C) oblique layers. The order of convergence is determined by using a linear fitting in a least squares sense. Table 2 Errors and experimental order of convergence in different situations: (A) horizontal layers, (B) vertical layers and (C) oblique layers. Here, Eε denotes the error in pressure in the discrete L2 -norm and p is an estimate of the convergence error.

ε −1

(A) 104 · Eε (m)

(B) 104 · Eε (m)

(C) 104 · Eε (m)

5 10 15 20 25 30

0.91 (−) 0.46 (0.99) 0.31 (0.92) 0.23 (1.06) 0.19 (0.88) 0.16 (1.10)

2.61 (−) 1.30 (1.00) 0.86 (1.03) 0.63 (1.05) 0.50 (1.04) 0.42 (0.98)

1.60 (−) 0.80 (1.00) 0.55 (0.91) 0.44 (0.80) 0.35 (0.94) 0.29 (1.13)

Thus, the angle defining the layer-structured direction is γ = −π /3 and, up to the associated change of coordinates, we have defined



K2 (y) =

10, 1,

if 0.0 ≤ y < 0.5, if 0.5 ≤ y < 1.0.

In the corresponding homogenized setting, this problem leads to the following homogenized coefficients:

κ1⋆ =

11

,

κ2⋆ =

20

.

2 11 Figs. 9–11 show the distribution of the pressure (isobaric lines) and saturation in the dam domain for different values of the period, namely ε = 1/10, 1/30 and in the homogenized case. Computing the isotropic homogenized solution leads us to consider the following permeability matrix and gravity term, with g = g (0, −1):

A0 ≃



3.11 1.59

1.59 , 4.21



B 0 (Θ , Θ ) · g ≃ −Θ

1.59 . 4.21





Numerical results highlight the fact that the ε -dependent pressure converges to the homogenized isotropic pressure, as observed in Figs. 9–11. This is also discussed in the next subsection, dealing with error analysis. 4.2. Numerical order of convergence in ε In Table 2 the convergence of the pressure of the ε -dependent problem to the one obtained in the homogenized problem can be numerically illustrated by means of the error Eε := ∥pε − p0 ∥L2 (Ω ) . In this way, the order of convergence can be numerically estimated. Assuming that Eε = O (ε m ), we show that the convergence is linear, i.e. m = 1, as illustrated by Fig. 12. Acknowledgments This paper has been partially funded by MCINN (Project MTM2010-21135-C02-01) and Xunta de Galicia (Ayuda CN2011/004, cofunded with FEDER funds).

S. Martin, C. Vázquez / Nonlinear Analysis: Real World Applications 14 (2013) 2133–2151

2151

References [1] H. Brezis, D. Kinderlehrer, G. Stampacchia, Sur une nouvelle formulation du problème de l’écoulement à travers une digue, C. R. Acad. Sci. Paris Sér. A–B 287 (9) (1978) A711–A714. [2] H.W. Alt, Strömungen durch inhomogene poröse Medien mit freiem Rand, J. Reine Angew. Math. 305 (1979) 89–115. [3] J.-F. Rodrigues, Some remarks on the homogenization of the dam problem, Manuscripta Math. 46 (1–3) (1984) 65–82. [4] D. Cioranescu, A. Damlamian, G. Griso, Periodic unfolding and homogenization, C. R. Math. Acad. Sci. Paris 335 (1) (2002) 99–104. [5] A. Bermúdez, J. Durany, Numerical solution of steady-state flow through a porous dam, Comput. Methods Appl. Mech. Engrg. 68 (1) (1988) 55–65. [6] C. Baiocchi, Free boundary problems in fluid flow through porous media and variational inequalities, in: Free Boundary Problems, Vol. I (Pavia, 1979), in: Ist. Naz. Alta Mat. Francesco Severi, Rome, 1980, pp. 175–191. [7] C. Baiocchi, V. Comincioli, E. Magenes, G. Pozzi, Free boundary problems in the theory of fluid flow through porous media: existence and uniqueness theorems, Ann. Mat. Pura Appl. (4) 97 (1973) 1–82. [8] C. Baiocchi, A. Friedman, A filtration problem in a porous medium with variable permeability, Ann. Mat. Pura Appl. (4) 114 (1977) 377–393. [9] V. Benci, On a filtration problem through a porous medium, Ann. Mat. Pura Appl. (4) 100 (1974) 191–209. [10] H. Darcy, Les fontaines publiques de la ville de Dijon, first ed., Dalmont, Paris, 1856. [11] J. Carrillo, M. Chipot, On the dam problem, J. Differential Equations 45 (2) (1982) 234–271. [12] H.W. Alt, G. Gilardi, The behavior of the free boundary for the dam problem, Ann. Sc. Norm. Super. Pisa Cl. Sci. (4) 9 (4) (1982) 571–626. [13] H.W. Alt, The fluid flow through porous media. Regularity of the free surface, Manuscripta Math. 21 (3) (1977) 255–272. [14] R. Stavre, B. Vernescu, Incompressible fluid flow through a nonhomogeneous and anisotropic dam, Nonlinear Anal. 9 (8) (1985) 799–810. [15] A. Friedman, S.Y. Huang, The inhomogeneous dam problem with discontinuous permeability, Ann. Sc. Norm. Super. Pisa Cl. Sci. (4) 14 (1) (1987) 49–77. [16] H.W. Alt, Numerical solution of steady-state porous flow free boundary problems, Numer. Math. 36 (1) (1980–1981) 73–98. [17] L.D. Marini, P. Pietra, Fixed-point algorithms for stationary flow in porous media, Comput. Methods Appl. Mech. Engrg. 56 (1) (1986) 17–45. [18] A. Bermúdez, J. Durany, Solving some free boundary problems by a method of characteristics, in: Free Boundary Problems: Theory and Applications, Vol. II (Irsee, 1987), in: Pitman Res. Notes Math. Ser., vol. 186, Longman Sci. Tech., Harlow, 1990, pp. 574–581. [19] A. Chakib, T. Ghemires, A. Nachaoui, A numerical study of filtration problem in inhomogeneous dam with discontinuous permeability, Appl. Numer. Math. 45 (2–3) (2003) 123–138. [20] A.R. Nazemi, M.H. Farahi, M. Zamirian, Filtration problem in inhomogeneous dam by using embedded method, J. Appl. Math. Comput. 28 (1–2) (2008) 313–332. [21] J. Bear, Hydraulics of Groundwater, McGraw-Hill, New York, 1979. [22] J. Carrillo, M. Chipot, The dam problem with leaky boundary conditions, Appl. Math. Optim. 28 (1) (1993) 57–85. [23] M. Chipot, A. Lyaghfouri, The dam problem for linear Darcy’s law and nonlinear leaky boundary conditions, Adv. Differential Equations 3 (1) (1998) 1–50. [24] A. Lyaghfouri, A free boundary problem for a fluid flow in a heterogeneous porous medium, Ann. Univ. Ferrara Sez. VII (NS) 49 (2003) 209–262. [25] S. Challal, A. Lyaghfouri, A filtration problem through a heterogeneous porous medium, Interfaces Free Bound. 6 (1) (2004) 55–79. [26] S. Challal, A. Lyaghfouri, The heterogeneous dam problem with leaky boundary condition, Commun. Pure Appl. Anal. 10 (1) (2011) 93–125. [27] D. Cioranescu, A. Damlamian, G. Griso, The periodic unfolding method in homogenization, SIAM J. Math. Anal. 40 (4) (2008) 1585–1620. [28] G. Nguetseng, A general convergence result for a functional related to the theory of homogenization, SIAM J. Math. Anal. 20 (3) (1989) 608–623. [29] G. Allaire, Homogenization and two-scale convergence, SIAM J. Math. Anal. 23 (6) (1992) 1482–1518. [30] D. Lukkassen, G. Nguetseng, P. Wall, Two-scale convergence, Int. J. Pure Appl. Math. 2 (1) (2002) 35–86. [31] M. Chipot, Variational Inequalities and Flow in Porous Media, in: Applied Mathematical Sciences, vol. 52, Springer-Verlag, New York, 1984. [32] M. Codegone, J.-F. Rodrigues, On the homogenization of the rectangular dam problem, Rend. Semin. Mat. Univ. Politec. Torino 39 (2) (1981) 125–136. 1982. [33] J. Nečas, Introduction to the Theory of Nonlinear Elliptic Equations, A Wiley-Interscience Publication, John Wiley & Sons Ltd., Chichester, 1986. Reprint of the 1983 edition. [34] A. Bermúdez, J. Durany, Numerical solution of cavitation problems in lubrication, Comput. Methods Appl. Mech. Engrg. 75 (1–3) (1989) 457–466. [35] G. Bayada, M. Chambat, C. Vázquez, Characteristics method for the formulation and computation of a free boundary cavitation problem, J. Comput. Appl. Math. 98 (2) (1998) 191–212. [36] G. Bayada, S. Martin, C. Vázquez, Homogenization of a nonlocal elastohydrodynamic lubrication problem: a new free boundary model, Math. Models Methods Appl. Sci. 15 (12) (2005) 1923–1956. [37] H. Brezis, Analyse fonctionnelle, in: Collection Mathématiques Appliquées pour la Maîtrise [Collection of Applied Mathematics for the Master’s Degree], in: Théorie et Applications. [Theory and applications], Masson, Paris, 1983. [38] A. Bermúdez, C. Moreno, Duality methods for solving variational inequalities, Comput. Math. Appl. 7 (1) (1981) 43–58.