Bingham flow in porous media with obstacles of different size

TL;DR

This paper develops a homogenization framework for Bingham fluid flow in porous media with obstacles of varying sizes, leading to a nonlinear Darcy law describing the effective flow behavior.

Contribution

It introduces a general compactness result for perforated domains with different-sized obstacles and applies it to derive a nonlinear Darcy law for Bingham fluids.

Findings

Homogenization of Bingham flow in media with mixed obstacle sizes

Derivation of a nonlinear Darcy law as the effective model

Establishment of a compactness result using unfolding operators

Abstract

By using the unfolding operators for periodic homogenization, we give a general compactness result for a class of functions defined on bounded domains presenting perforations of two different size. Then we apply this result to the homogenization of the flow of a Bingham fluid in a porous medium with solid obstacles of different size. Next we give the interpretation of the limit problem in term of a non linear Darcy law.

Equations

$$Y^{\ast }=Y\setminus Y_s,Z^{\ast }=Z\setminus Z_s,$$

$$\underset{\epsilon \to 0}{\lim }\delta (\epsilon )=0.$$

$$\chi_{Y^{\ast}}\left(y\right)=\left\{\begin{array}[c]{l}1,\ \text{in }Y^{\ast},\\ 0,\ \text{in }Y\setminus Y^{\ast}\end{array}\right.,\ \ \ \chi_{Z^{\ast}}\left(y\right)=\left\{\begin{array}[c]{l}1,\ \text{in }Z^{\ast},\\ 0,\ \text{in }Z\setminus Z^{\ast}.\end{array}\right.$$

$${\Omega }_{\epsilon \delta }=\left\{x\in \Omega :{\chi }_{Y^{\ast }}\left(\frac{x}{\epsilon }\right){\chi }_{Z^{\ast }}\left(\frac{x}{\epsilon \delta (\epsilon )}\right)=1\right\}.$$

$$x=\varepsilon\Bigl{[}\frac{x}{\varepsilon}\Bigr{]}_{Y}+\varepsilon\Bigl{\{}\frac{x}{\varepsilon}\Bigr{\}}_{Y},$$

$${\hat{K}}_{\epsilon }=\{k\in {\mathbb{Z}}^N|\epsilon (k+Y)\subset \Omega \},{\hat{\Omega }}_{\epsilon }=\mathrm{int}\bigcup _{k\in {\hat{K}}_{\epsilon }}\left(\epsilon (k+\overline{Y})\right),{\Lambda }_{\epsilon }=\Omega \setminus {\hat{\Omega }}_{\epsilon },$$

$$T_{\varepsilon}(\varphi)(x,y)=\left\{\begin{array}[c]{ll}\varphi\left(\varepsilon\left[\dfrac{x}{\varepsilon}\right]_{Y}+\varepsilon y\right)&\text{for a.e. }(x,y)\in\widehat{\Omega}_{\varepsilon}\times Y,\\[8.61108pt] 0&\text{for a.e. }(x,y)\in\Lambda_{\varepsilon}\times Y.\end{array}\right.$$

$$T_{\epsilon }({\phi }_{{|}_{\Omega }}^{\epsilon })\to \phi \text{ strongly in }{L}^2(\Omega \times Y);$$

$$T_{\epsilon }({\phi }_{\epsilon })\to \phi \text{ strongly in }{L}^2(\Omega \times Y).$$

$$T_{\epsilon }({\phi }_{\epsilon })\to \phi \text{ weakly in }{L}^2(\Omega \times Y).$$

$${\phi }_{\epsilon }\to {\mathcal{M}}_Y(\phi )\text{ weakly in }{L}^2(\Omega ),$$

$${\mathcal{M}}_Y(\phi )=\frac{1}{\left|Y\right|}{\int }_Y\phi \left(x,y\right)dy\text{ a.e. for }x\in \Omega .$$

$${\nabla }_y(T_{\epsilon }(\phi ))=\epsilon T_{\epsilon }(\nabla \phi )\text{ a.e. }(x,y)\in \Omega \times Y.$$

$$T_{\delta}(v)(x,y,z)=v\left(x,\delta\Bigl{[}\frac{y}{\delta}\Bigr{]}_{Z}+\delta z\right)\quad\text{ for }(x,y,z)\in\Omega\times Y\times Z,$$

$$T_{\delta}(T_{\varepsilon}(\varPhi))(x,y,z)=\varPhi\Bigl{(}\varepsilon\Bigl{[}\frac{x}{\varepsilon}\Bigr{]}_{Y}+\varepsilon\delta\Bigl{[}\frac{y}{\delta}\Bigr{]}_{Z}+\varepsilon\delta z\Bigr{)}.$$

$$\epsilon \delta T_{\delta }(T_{\epsilon }({\nabla }_x\Phi ))=T_{\delta }(T_{\epsilon }(\epsilon \delta {\nabla }_x\Phi ))=T_{\delta }(\delta T_{\epsilon }(\epsilon {\nabla }_x{\Phi }_{\epsilon \delta }))=T_{\delta }(\delta ({\nabla }_y{T}_{\epsilon }(\Phi ))={\nabla }_z{T}_{\delta }(T_{\epsilon }(\Phi ))$$

$$\epsilon \delta {\Vert \nabla {\Phi }_{\epsilon \delta }\Vert }_{(L^2\left(\Omega \right){)}^n}\le C.$$

$$T_{\delta }\left(T_{\epsilon }({\Phi }_{\epsilon \delta })\right)$$

$$T_{\epsilon }({\Phi }_{\epsilon \delta })$$

$${\Phi }_{\epsilon \delta }$$

$$\epsilon \delta T_{\delta }\left(T_{\epsilon }({\nabla }_x{\Phi }_{\epsilon \delta })\right)$$

$$\epsilon \delta T_{\delta }(T_{\epsilon }({\nabla }_x{\Phi }_{\epsilon \delta }))={\nabla }_z{T}_{\delta }(T_{\epsilon }({\Phi }_{\epsilon \delta }))$$

$${\int }_{\Omega }\Phi dx\sim \frac{1}{\left|Y\right|\left|Z\right|}{\int }_{\Omega \times Y\times Z}{T}_{\delta }(T_{\epsilon }\left(\Phi \right))dxdydz,\forall \phi \in L^1\left(\Omega \right),$$

$${\sigma }_{ij}=-p_{\epsilon \delta }{\delta }_{ij}+g\epsilon \delta \frac{D_{ij}(u_{\epsilon \delta })}{(D_{II}(u_{\epsilon \delta }){)}^{\frac{1}{2}}}+2\mu {\epsilon }^2{\delta }^2{D}_{ij}(u_{\epsilon \delta }),$$

$$D_{ij}(u_{\epsilon \delta })$$

$$D_{II}(u_{\epsilon \delta })$$

$${\sigma }_{ij}^D$$

$${\sigma }_{II}$$

$$\left\{\begin{array}[c]{rcl}(\sigma_{II})^{\frac{1}{2}}<g\varepsilon\delta&\Leftrightarrow&D_{ij}(u_{\varepsilon\delta})=0\\ (\sigma_{II})^{\frac{1}{2}}\geq g\varepsilon\delta&\Leftrightarrow&D_{ij}(u_{\varepsilon\delta})=\dfrac{1}{2\mu\varepsilon^{2}\delta^{2}}\left(1-\dfrac{g\varepsilon\delta}{(\sigma_{II}^{\varepsilon})^{\frac{1}{2}}}\right){\sigma}_{ij}^{D}.\end{array}\right.$$

$$\mathrm{div}{u}_{\epsilon \delta }=0\text{in}{\Omega }_{\epsilon \delta }.$$

$$\left\{\begin{array}[c]{l}a_{\varepsilon\delta}\left(u_{\varepsilon\delta},v-u_{\varepsilon\delta}\right)+j_{\varepsilon\delta}\left(v\right)-j_{\varepsilon\delta}\left(u_{\varepsilon\delta}\right)\geq\left(f,v-u_{\varepsilon\delta}\right)_{\Omega_{\varepsilon\delta}},\ \forall v\in V\left(\Omega_{\varepsilon\delta}\right)\\ u_{\varepsilon\delta}\in V\left(\Omega_{\varepsilon\delta}\right),\end{array}\right.$$

$$a_{\epsilon \delta }\left(u,v\right)$$

$$V\left({\Omega }_{\epsilon \delta }\right)$$

$$\left\{\begin{array}[c]{l}a_{\varepsilon\delta}\left(u_{\varepsilon\delta},v-u_{\varepsilon\delta}\right)+j_{\varepsilon\delta}\left(v\right)-j_{\varepsilon\delta}\left(u_{\varepsilon\delta}\right)\geq\left(f,v-u_{\varepsilon\delta}\right)_{\Omega_{\varepsilon\delta}}-\langle\nabla p_{\varepsilon\delta},v-u_{\varepsilon\delta}\rangle_{\Omega_{\varepsilon\delta}},\ \forall v\in\left(H_{0}^{1}\left(\Omega_{\varepsilon\delta}\right)\right)^{n}\\ u_{\varepsilon\delta}\in V\left(\Omega_{\varepsilon\delta}\right),\ p_{\varepsilon\delta}\in L_{0}^{2}\left(\Omega_{\varepsilon\delta}\right)\end{array}\right.$$

Full text

Keywords: homogenization, unfolding operators, Bingham fluid, porous media

MSC: 35B27, 76M50, 76S05

1 Introduction

In this paper we study the homogenization problem for a Bingham flow in a porous medium with solid obstacles of different size. The aim of our paper is twofold: we first define the unfolding operators for periodic homogenization in a domain which presents periodically distributed perforations of two different size and we give corresponding compactness results. Then we illustrate these results with an application to the homogenization of a Bingham flow in a porous medium with solid obstacles of different size.

In order to define the appropriate unfolding operators and to get the compactness results, we follow the ideas introduced by D. Cioranescu, A. Damlamian and G. Griso in [8] and [17] for the case of functions with one scale of periodicity and developed later by A. Damlamian, N. Meunier, J. Van Schaftingen in [10] and [15] for the case of functions with more than one periodicity scales. Nevertheless, our result is different from the ones presented in the previous cited papers, due to the presence of the perforations at the two different scales. The case, different from the one presented here, corresponding to the unfolding operators for a doubly periodic domain presenting perforations at the very small scale only, was recently addressed by Bunoiu, and Donato in [5].

More precisely, our domain contains small perforations of size $\varepsilon$ periodically distributed with period $\varepsilon$ and very small perforations of size $\varepsilon\delta(\varepsilon)$ periodically distributed with periodicity $\varepsilon\delta(\varepsilon)$. Here $\varepsilon$ and $\delta(\varepsilon)$ are real positive parameters smaller than one with $\delta(\varepsilon)$ tending to zero when $\varepsilon$ tends to zero. Such a geometry modelizes, for example, a porous medium in which the perforations correspond to solid impervious obstacles.

In the fluid part of this porous medium we consider the stationary flow of the Bingham fluid, under the action of external forces. The Bingham fluid is an incompressible fluid which has a non linear constitutive law; so it is a non-Newtonian fluid. This fluid moves like a rigid body when a certain function of the stress tensor is below a given threshold. Beyond this threshold the fluid flows, obeying a non linear constitutive law. As example of such fluids we can mention some paints, the mud which can be used for the oil extraction and the volcanic lava. Bingham flow in other contexts is studied by Bunoiu, and Kesavan in[6]. For a presentation of the different types of non-Newtonian fluids we refer the reader to [9].

The mathematical model of the Bingham flow in a bounded domain was introduced in [11] by G. Duvaut and J. L. Lions. The existence of the velocity and of the pressure for this model was proved in the case of a bi-dimensional and of a three-dimensional domain.

The homogenization problem in a classical porous medium, with obstacles of size $\varepsilon$ and $\varepsilon$-periodically distributed, was first studied in [14] by J. L. Lions and E. Sanchez-Palencia. The authors did the asymptotic study of the problem by using a multiscale method, involving a “macroscopic” variable $x$ and a “microscopic” variable $y=\dfrac{x}{\varepsilon}$, associated to the relative dimension of the pores. The study is based on a multiscale “ansatz”, which allows to obtain to the limit a non linear Darcy law. There is no convergence result proved.

The rigorous justification for the convergence of the homogenization process of the results presented in [14] is given by A. Bourgeat and A. Mikelic in [3]. In order to do it, the authors used monotonicity methods coupled with the two-scale convergence method introduced by G. Nguetseng in [16] and further developed by G. Allaire in a series of papers, as for example [1]. The limit problem announced in [14] was obtained, by letting the small parameter $\varepsilon$ tend to zero in the initial problem. The unfolding method for periodic homogenization, introduced by D. Cioranescu, A. Damlamian and G. Griso in [8] was used by R. Bunoiu, G. Cardone, C. Perugia in [4] in order to obtain the limit problem. This method presents the advange of transforming in an easy manner the initial problem, stated in a domain dependent on $\varepsilon$, in a problem stated in a domain independent of $\varepsilon$. The passage to the limit when $\varepsilon$ tends to zero is then simple thanks to the compactness results, and this for the non linear terms too.

Our paper is organized as follows. In section 2 we define the double perforated domain, the unfolding operators adapted to it and then we give a compactness result. In section 3 we describe the problem of the Bingham flow and we give the a priori estimates for the velocity and the pressure of the flow. Following the ideas in R. Bunoiu, J. Saint Jean Paulin [7], we construct the extension of the pressure to the whole domain, namely in the perforations too. In section 4 we state the main result of the paper, which consists in getting the limit problem. It is obtained in two steps: we first apply the unfolding operators for periodic homogenization defined in section 2 to the variational formulation of the problem which describes the Bingham flow in our porous medium. Then we pass to the limit when the small parameter $\varepsilon$ tends to zero. In section 5 we give the interpretation of the limit problem in term of a non linear Darcy law and we compare it with the classical linear Darcy law.

2 Unfolding operators and compactness results

Let $\Omega$ be a bounded open domain in $\mathbb{R}^{n}$ with Lipschitz continuous boundary $\partial\Omega,$ $n=2$ or $n=3$.

We consider two fixed reference cells $Y=\left]0,y_{1}^{0}\right[\times...\times\left]0,y_{n}^{0}\right[$ and $Z=\left]0,z_{1}^{0}\right[\times...\times\left]0,z_{n}^{0}\right[$ and two closed subsets $Y_{s}$ and $Z_{s}$ with non-empty interior and Lipschitz continuous boundaries, contained in $Y$ and $Z$ respectively. We define:

[TABLE]

and we give a simple example in Figure 1.

Let $\varepsilon$ be a positive parameter, smaller than one. For every $\varepsilon>0$, let $0<\delta(\varepsilon)<\varepsilon$ be such that

[TABLE]

We suppose that there exists an $\varepsilon$ such that the domain $\overline{Y}$ is exactly covered by a finite number of cells $\delta(\varepsilon)\overline{Z}$. Moreover, we suppose that $Y_{s}$ is exactly covered by a finite number of cells $\delta(\varepsilon)\overline{Z}$. This last hypothesis implies some restrictions for the geometry of $Y_{s}$. We deduce that there is no intersection between the domains $Y_{s}$ and $\delta(\varepsilon)Z_{s}$ in the cell $Y$, as one can see on an example in Figure 2. If we consider all the small parameters $\varepsilon/2^{N}$ (with $N$ natural number), the above assumptions are still true. We denote $Y_{f}$ the complement in $Y$ of the set ${Y_{s}}\bigcup\limits_{l\in{\mathbb{Z}}^{n}}({\delta(\varepsilon)(lz^{0}+{Z}_{s})}\cap Y^{\ast})$.

We multiply the perforated cell $Y$ (Figure 2) by $\varepsilon$ and we repeat it in the domain $\Omega$. For simplicity and without loosing any generality, one could even assume that $\overline{\Omega}$ is exactly covered by a finite number of cells $\varepsilon\overline{Y}$. The domain $\Omega_{\varepsilon\delta}$ is the one obtained by taking out of $\Omega$ the translated of the domains $\varepsilon Y_{s}$ and $\varepsilon\delta(\varepsilon)Z_{s}$. Let us notice that there is no intersection between the solid obstacles $\varepsilon Y_{s}$ and $\varepsilon\delta(\varepsilon)Z_{s}$ in $\Omega_{\varepsilon\delta}$, because there is no intersection between the solid obstacles $Y_{s}$ and $\delta(\varepsilon)Z_{s}$ in the cell $Y$. The domain $\Omega_{\varepsilon\delta}$ is connected, but the union of solid obstacles is not connected (see an example in Figure 2).

Let $\chi_{Y^{\ast}}$ and $\chi_{Z^{\ast}}$ be the characteristic functions of the domains $Y^{\ast}$and $Z^{\ast}$, defined by:

[TABLE]

We extend the characteristic functions $\chi_{Y^{\ast}}$ (respectively $\chi_{Z^{\ast}}$) by periodicity, with period $y_{i}^{0}$ in $y_{i}$ and with period $z_{i}^{0}$ in $z_{i}$, for $i=1,..,n$. The domain $\Omega_{\varepsilon\delta}$, defined as above is described by:

[TABLE]

The domain $\Omega_{\varepsilon\delta}$ presents a structure with a double periodicity: there are small perforations of size $\varepsilon$ and very small perforations of size $\varepsilon\delta(\varepsilon)$. The boundary $\partial\Omega_{\varepsilon\delta}$ that is composed by two parts: the boundary of $\Omega$, denoted $\partial\Omega$ and the union of the boundaries of all the obstacles, denoted $\Gamma_{\varepsilon}.$

We follow the general idea of the unfolding method, namely we transform oscillating functions defined on the domain $\Omega$ into functions defined on the domain $\Omega\times Y\times Z$. In order to do this, we proceed in two steps: first we use the general theory of the unfolding homogenization in order to make the transformation from the domain $\Omega$ to the domain $\Omega\times Y$. In order to do this, we use the unfolding operator introduced in [8] for the scale $\varepsilon$. Next we define a second unfolding operator, for the scale $\delta(\varepsilon)$, which allows us to transform oscillating functions defined on $\Omega\times Y$ into functions defined on $\Omega\times Y\times Z$. In order to do this, we will follow in addition the ideas of A. Damlamain N. Meunier, J. Van Schaftingen in [10] and [15].

For the first step, the idea is to transform oscillating functions defined on the domain $\Omega$ into functions defined on the domain $\Omega\times Y$, in order to isolate the oscillations in the second variable. This transformation, together with a priori estimates, allows us to use compactness results and then to get the limits of our oscillating sequences. We start by recalling the results as far as the unfolding operator for the scale $\varepsilon$ is concerned.

We know that every real number $a$ can be written as the sum between his integer part $[a]$ and his fractionary part $\left\{a\right\}$ which belongs to the interval $\left[0,1\right)$.

For $x=(x_{1},...,x_{n})\in\mathbb{R}^{n}$, we apply a similar decomposition to every real number $\dfrac{x_{i}}{\varepsilon}$ for $i=1,...,n$ and we get

[TABLE]

where \Bigl{[}\dfrac{x}{\varepsilon}\Bigr{]}_{Y}\in\mathbb{Z}^{n} and \Bigl{\{}\dfrac{x}{\varepsilon}\Bigr{\}}_{Y}\,\in Y.

We define now

[TABLE]

and we notice that the set $\widehat{\Omega}_{\varepsilon}$ is the interior of the largest union of $\varepsilon(k+\overline{Y})$ cells included in $\Omega$.

Definition 1

For any Lebesgue measurable function $\varphi$ on $\Omega$, we define the periodic unfolding operator by the formula

[TABLE]

According to [8], this operator has the following properties:

1. p${}_{1})$

$T_{\varepsilon}$ is linear and continuous from $L^{2}(\Omega)$ to $L^{2}(\Omega\times Y)$; 2. p${}_{2})$

$T_{\varepsilon}(\varphi\,\phi)=T_{\varepsilon}(\varphi)T_{\varepsilon}(\phi),$ $\forall\varphi,\phi\in L^{2}(\Omega);$ 3. p${}_{3})$

If $\varphi\in L^{2}(Y)$ is a $Y$-periodic function and \varphi^{\varepsilon}(x)=\varphi\Bigl{(}\displaystyle\frac{x}{\varepsilon}\Bigr{)},\,x\,\in\,\mathbb{R}^{N} then

[TABLE] 4. p${}_{4})$

If $\varphi_{\varepsilon}\in L^{2}(\Omega)$ and $\varphi_{\varepsilon}\rightarrow\varphi$ strongly in $L^{2}(\Omega),$ then

[TABLE]

Moreover, the following results hold (see Proposition 2.9 (iii) in [8]):

Proposition 2

Let $\{\varphi_{\varepsilon}\}_{\varepsilon}$ be a bounded sequence in $L^{2}(\Omega)$ such that

[TABLE]

Then

[TABLE]

where the mean value operator* $\mathcal{M}_{Y}$ is defined by*

[TABLE]

Let us moreover observe that for a function $\varphi\,\in H^{1}(\Omega)$, one has

[TABLE]

We define now the second unfolding operator, at the scale $\delta(\varepsilon)$ (denoted in the sequel by $\delta$).

Definition 3

Let $v\in L^{p}(\Omega\times Y)$. Then the unfolding operator $T_{\delta}$ is defined by

[TABLE]

where $x$ plays the role of a parameter.

Let now $\varPhi$ be a function belonging to the space $H^{1}_{0}(\Omega)$. Then, accordind to Definitions 1 and 3 we have

[TABLE]

Moreover, the following equality holds true:

[TABLE]

and we have the convergence results:

Proposition 4

Let $\left\{\varPhi_{\varepsilon\delta}\right\}_{\varepsilon\delta}$ be a sequence in $H_{0}^{1}(\Omega)$ bounded in $L^{2}(\Omega)$. Let us assume that

[TABLE]

Then there exists $\widehat{\varPhi}$ in $L^{2}(\Omega\times Y;H^{1}_{\#}\left(Z\right))$ such that, up to a subsequence still denoted by $\varepsilon\delta$ we have

[TABLE]

where $\displaystyle H^{1}_{\#}\left(Z\right)=\left\{\phi\in H^{1}(Z),\phi\hbox{ is Z-periodic}\right\}.$

Proof. The sequence $\left\{\varPhi_{\varepsilon\delta}\right\}_{\varepsilon\delta}$ being bounded in $H_{0}^{1}(\Omega)$, then $\left\{T_{\delta}\left(T_{\varepsilon}(\varPhi_{\varepsilon\delta})\right)\right\}_{\varepsilon\delta}$ is bounded in the space $L^{2}(\Omega\times Y\times Z)$. Clasical compactness results imply the existence of a function $\widehat{\varPhi}$ in $L^{2}(\Omega\times Y\times Z)$ such that the first weak convergence holds true. By using Proposition 2 and its analoguous at the scale $\delta$, we obtain the second and the third weak convergences.

The last weak converge is a consequence of the equality

[TABLE]

and of the first weak convergence, for the sequence $T_{\delta}(T_{\varepsilon}(\nabla_{x}\varPhi_{\varepsilon\delta})).$ The fact that the limit $\widehat{\varPhi}$ actually belongs to the space $L^{2}(\Omega\times Y;H^{1}_{\#}\left(Z\right))$ is due to the application of a result from [8] to the unfolding operator at the scale $\delta$. $\blacksquare$

Remark 5

If $\left\{\varPhi_{\varepsilon\delta}\right\}_{\varepsilon\delta}$ is a sequence in $H_{0}^{1}(\Omega_{\varepsilon\delta})$ those extension by zero to the whole of $\Omega$ satisfy the hypothesis of Proposition 4, then all the results still hold true, with $Y$ and $Z$ replaced by $Y^{\ast}$ and $Z^{\ast}$ respectively in the description of all the function spaces and in the integrals. Indeed, the sequence $\left\{\varPhi_{\varepsilon\delta}\right\}_{\varepsilon\delta}$ vanishes on $\Omega\times Y_{s}\times Z_{s}$ and this property is preserved by passing to the limit.

To end this section, we recall one of the key points of the use of the unfolding method for periodic homogenization: the fact that the integrals over the domain $\Omega$ can be replaced by integrals over the domain $\Omega\times Y\times Z$, by using the relation below

[TABLE]

which is true for $\varepsilon$ and $\delta$ sufficiently small.

3 Statement of the problem and preliminary results

Our aim now is to apply the results from Section 2 to the homogenization of a problem stated in a domain $\Omega_{\varepsilon\delta}$ defined as before. The problem we address is the flow of a Bingham fluid in a porous medium with obstacles of different size. Indeed, such porous media can be modelized by the domain $\Omega_{\varepsilon\delta}$, corresponding to the part where the fluid flows. The perforations correspond to solid impervious obstacles. If $u_{\varepsilon\delta}$ and $p_{\varepsilon\delta}$ are the velocity and pressure respectively for a Bingham fluid, then its stress tensor is defined by

[TABLE]

where $\delta_{ij}$ is the Kronecker symbol, $g$ and $\mu$ are real positive constants. The constant $g$ represents the yield stress of the fluid and the constant $\mu$ is the viscosity. Relation (3.1) represents the constitutive law of the Bingham fluid.

We define the entries of the strain tensor, denoted $D(u_{\varepsilon\delta})$, by

[TABLE]

Let us note that the constitutive law (3.1) is valid only if $D_{II}(u_{\varepsilon\delta})\neq 0.$ In [11] it is shown that this constitutive law is equivalent with the following one:

[TABLE]

We see that this is a threshold law: as long as the shear stress is below $g\varepsilon\delta$, the fluid behaves as a rigid solid. When the value of the shear stress exceeds $g\varepsilon\delta$, the fluid flows and obeys a non linear law.

Moreover, the fluid is incompressible, which means that its velocity is divergence free

[TABLE]

In [11] it is shown that the velocity $u_{\varepsilon\delta}$ satisfies the following variational inequality when we apply to the porous media an external force denoted by $f$ and belonging to $\left(L^{2}\left(\Omega\right)\right)^{n}$:

[TABLE]

where

[TABLE]

If $f\in\left(L^{2}\left(\Omega\right)\right)^{n}$, we know from [11] that for $n=2$ or $3$ and every fixed $\varepsilon$ and $\delta$ there exists a unique $u_{\varepsilon\delta}\in V\left(\Omega_{\varepsilon\delta}\right)$ solution of problem (3.2) and that if $p_{\varepsilon\delta}$ is the pressure of the fluid in $\Omega_{\varepsilon\delta}$, then the problem (3.2) is equivalent to the following one:

[TABLE]

Here $L_{0}^{2}\left(\Omega_{\varepsilon\delta}\right)$ denotes the space of functions belonging to $L^{2}\left(\Omega_{\varepsilon\delta}\right)$ and of mean value zero. For an open set $D$, the brackets $\langle\cdot,\cdot\rangle_{D}$ denote the duality product between the spaces $H^{-1}(D)^{n}$ and $H^{1}_{0}(D)^{n}$, where $H^{-1}(D)^{n}$ denotes the dual of $H^{1}_{0}(D)^{n}$.

Our aim now is to pass to the limit as $\varepsilon\rightarrow 0$ and $\delta\rightarrow 0$ in problem (3.3). In order to do this, we first need to get a priori estimates for the velocity $u_{\varepsilon\delta}$ and the pressure $p_{\varepsilon\delta}$. An important role is played by the value of the constant in Poincaré’s inequality, with reads:

Proposition 6

Let $v$ be a function in $\left(H_{0}^{1}\left(\Omega_{\varepsilon\delta}\right)\right)^{n}$. Then we have the following inequality:

[TABLE]

Proof. We prove this result by using a crucial result of Tartar (see [17]), that we generalize here to the case of a domain with two scales of periodicity. The idea is to derive Poincaré’s in the whole domain by succesively using the $\varepsilon$-periodicity and $\delta$-periodicity of the domain respectively and by applying the classical Poincaré inequality in the cell $Z^{\ast}$.

More precisely, due to the $\varepsilon$-periodicity, it is clear that we have:

[TABLE]

and

[TABLE]

In this above sum there are $N_{\varepsilon}$ terms and by construction $\displaystyle N_{\varepsilon}=\frac{|\Omega|}{|\varepsilon Y|}\approx\varepsilon^{-n}\frac{|\Omega|}{|Y|}.$

Therefore, in order to obtain the Poincaré inequality in the whole domain $\Omega_{\varepsilon\delta}$ it is enough to know it in an arbitrary cell $\varepsilon(ky^{0}+Y^{\ast})$ and then to sum over $k$.

For $v$ a function in $\left(H_{0}^{1}\left(\Omega_{\varepsilon\delta}\right)\right)^{n}$ we define

[TABLE]

This function is defined on $Y^{\ast}$, it belongs to the space $H^{1}(Y^{\ast})$ and $v_{\varepsilon,k^{{}^{\prime}}}(y)=0$ on $\partial Y_{s}.$

Moreover, due to the equalities

[TABLE]

and

[TABLE]

it is now enough to know the Poincaré inequality in the domain $Y^{\ast}$ in order to get the result.

Due to the $\delta$-periodicity and to the hypothesis on the geometry of our domain, we have

[TABLE]

and

[TABLE]

In this above sum there are $N_{\delta}$ terms and by construction $\displaystyle N_{\delta}=\frac{|Y^{\ast}|}{|\delta Z|}=\delta^{-n}\frac{|Y^{\ast}|}{|Z|}.$

Therefore, in order to obtain the Poincaré inequality in the domain $Y^{\ast}$ it is enough to know it in an arbitrary cell $\delta(lz^{0}+Z^{\ast})$ and then to sum over $l$. By using an argument as above it is actually enough to know the Poincaré inequality in the domain $Z^{\ast}.$ We define

[TABLE]

In $Z^{\ast}$ we know the classical Poincaré inequality:

[TABLE]

We point out that the constant $C(Z^{*})$ is independent on $\varepsilon$ and on $\delta$. This implies

[TABLE]

[TABLE]

By summing now over $l$ and then by repeating the same argument at the scale $\varepsilon$ and summing over $k$ we obtain the desired result. $\blacksquare$

Proposition 7

The solution $(u_{\varepsilon\delta},p_{\varepsilon\delta})$ of problem (3.3) satisfies the following a priori estimates:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Proof. Setting $v=2u_{\varepsilon\delta}$ and $v=0$ successively in (3.2) and using the Poincaré inequality, we find the first two estimates, for the velocity.

Let $v_{\varepsilon\delta}\in\left(H_{0}^{1}\left(\Omega_{\varepsilon\delta}\right)\right)^{n}$. Setting $v=v_{\varepsilon\delta}+u_{\varepsilon\delta}$ in (3.3) and using estimates on the velocity, we obtain the first estimate for the pressure and then we deduce the second one, by using a rescaled Nečas inequality. $\blacksquare$

Now we extend the velocity $u_{\varepsilon\delta}$ by zero to $\Omega\diagdown\Omega_{\varepsilon\delta}$, denote the extension by the same symbol and we have the following estimates:

[TABLE]

Moreover, we remark that $\operatorname{div}u_{\varepsilon\delta}=0$ in $\Omega$.

In order to define the extension of the pressure to the whole domain $\Omega$, we generalize here the results from R. Bunoiu, J. Saint Jean Paulin [7], which followed the classical idea of L. Tartar [17]. We first construct a restriction operator $S_{\varepsilon\delta}$ from $(H_{0}^{1}\left(\Omega\right))^{n}$ to $(H_{0}^{1}\left(\Omega_{\varepsilon\delta}\right))^{n}$ and using this operator we then define an extension for the pressure to the whole domain $\Omega$.

We define the spaces $H_{s}^{1}\left(Y^{\ast}\right)$ and $H_{s}^{1}\left(Y_{f}\right)$ by

[TABLE]

and

[TABLE]

where the domain $Y_{f}$ is defined in section 2 and an example is given in Figure 2. We denote $\Gamma$ the union of the boundaries of all the obstacles contained in $Y$.

Now we first construct a restriction operator $R$ from the space $(H^{1}(Y))^{n}$ into the space $(H_{s}^{1}\left(Y^{\ast}\right))^{n}$ and next we construct a second restriction operator $W_{\delta}$ from the space $(H_{s}^{1}\left(Y^{\ast}\right))^{n}$ into the space $(H_{s}^{1}\left(Y_{f}\right))^{n}$. Using the operators $R$ and $W_{\delta}$, we then construct the operator

[TABLE]

and finally we define $S_{\varepsilon\delta}$ by applying $S_{\delta}$ to each period $\varepsilon Y$ of $\Omega$. So we construct $S_{\varepsilon\delta}$ in three steps, corresponding to the three following lemmas.

Lemma 8

There exists a restriction operator

[TABLE]

such that for $v\in(H^{1}(Y))^{n}$ we have

$Rv=v$ if $v=0$ in $Y_{s};$ 2. 2.

$\operatorname{div}Rv=0$ in $Y^{\ast}$ if $\operatorname{div}v=0$ in $Y,$ 3. 3.

$\left\|Rv\right\|_{(H_{s}^{1}(Y^{\ast}))^{n}}\leq c\left\|v\right\|_{(H^{1}(Y))^{n}}.$

Lemma 9

There exists a restriction operator

[TABLE]

such that for $Rv\in(H_{s}^{1}(Y^{\ast}))^{n}$ we have

$W_{\delta}\left(Rv\right)=Rv$ if $Rv=0$ in $\bigcup\limits_{l\in{\mathbb{Z}}^{n}}({\delta(\varepsilon)(lz^{0}+{Z}_{s})}\bigcap Y^{*}$. 2. 2.

$\operatorname{div}W_{\delta}\left(Rv\right)=0$ in $Y_{f}$ if $\operatorname{div}Rv=0$ in $Y^{\ast},$ 3. 3.

$\delta\left\|\nabla W_{\delta}\left(Rv\right)\right\|_{\left(L^{2}(Y_{f})\right)^{n\times n}}+c\left\|W_{\delta}\left(Rv\right)\right\|_{(L^{2}(Y_{f}))^{n}}\leq c\left\|v\right\|_{(H^{1}(Y))^{n}}.$

Lemma 10

There exists a restriction operator

[TABLE]

such that

$S_{\varepsilon\delta}\left(v\right)=v$ in $\Omega_{\varepsilon\delta},$ $\forall v\in(H_{0}^{1}(\Omega_{\varepsilon\delta}))^{n},$ 2. 2.

$\operatorname{div}S_{\varepsilon\delta}v=0$ in $\Omega_{\varepsilon\delta}$ if $\operatorname{div}v=0$ in $\Omega,$ 3. 3.

$\displaystyle\varepsilon\delta\left\|\nabla S_{\varepsilon\delta}v\right\|_{\left(L^{2}(\Omega_{\varepsilon\delta})\right)^{n\times n}}+c\left\|S_{\varepsilon\delta}v\right\|_{(L^{2}(\Omega_{\varepsilon\delta}))^{n}}\leq c\left\|v\right\|_{(H^{1}(\Omega))^{n}}.$

Let now $v$ be a function in the space $H_{0}^{1}(\Omega).$ As $\nabla p_{\varepsilon\delta}\in H^{-1}\left(\Omega_{\varepsilon\delta}\right)$, we define the application $F^{\varepsilon}$ by

[TABLE]

where $S_{\varepsilon\delta}$ is the operator defined by Lemma 10. The following proposition defines us the extension $\widetilde{p}_{\varepsilon\delta}$ of the pressure $p_{\varepsilon\delta}$ to the whole $\Omega$. Moreover, it gives us a strong convergence result for this extension. Following the ideas of L.Tartar [17], we can prove

Proposition 11

Let $p_{\varepsilon\delta}$ be as in (3.3). Then, for each $\varepsilon$ and $\delta$ there exists an extension $\widetilde{p}_{\varepsilon\delta}$ of $p_{\varepsilon\delta}$ defined on $\Omega$ such that

[TABLE]

Moreover, up to a subsequence, we have

[TABLE]

The function $F^{\varepsilon}$ defined before and the pressure $\widetilde{p}_{\varepsilon\delta}$ are linked by the relation

[TABLE]

For every function $v$ that is the extension by zero to the whole $\Omega$ of a function in $H_{0}^{1}(\Omega_{\varepsilon\delta})^{n}$ we deduce:

[TABLE]

According to the extensions of the velocity and of the pressure, problem (3.3) can now be written as:

[TABLE]

for every $v$ that is the extension by zero to the whole $\Omega$ of a function in $H_{0}^{1}(\Omega_{\varepsilon\delta})^{n}$.

4 Convergence result

Now we state the main result of our paper, the convergence result for the variational inequality (3.4). In order to prove it, we apply the unfolding operators from section 2, together with the a priori estimates from Proposition 7 and the compactness results from Proposition 4.

Theorem 12

Let $u_{\varepsilon\delta}$ and $\widetilde{p}_{\varepsilon\delta}$ verify relation (3.4). Then there exist $\widehat{u}\left(x,y,z\right)\in L^{2}\left(\Omega\times Y^{\ast};(H_{per}^{1}\left(Z^{\ast}\right))^{n}\right)$ and $\widehat{p}\in L_{0}^{2}(\Omega)\cap H^{1}(\Omega)$ such that

[TABLE]

[TABLE]

and satisfy the limit problem

[TABLE]

for every $\Phi\in\mathcal{V},$ where

[TABLE]

[TABLE]

[TABLE]

The function $\widehat{u}$ satisfies the following conditions:

[TABLE]

Proof. Taking into account the a priori estimates from Proposition 7 and then using Proposition 4 and Remark 5, we have the following convergences for the velocity:

[TABLE]

According to [17], we have for the pressure the converegnce

[TABLE]

Using property $p_{4}$) of the unfolding operators we get:

[TABLE]

In order to prove relation (4.2), let us observe that $\operatorname{div}u_{\varepsilon\delta}=0$ implies $\varepsilon T_{\varepsilon}\left(\operatorname{div}u_{\varepsilon\delta}\right)=0.$ But

[TABLE]

and so

[TABLE]

which implies $\operatorname{div}_{z}T_{\delta}(T_{\varepsilon}\left(u_{\varepsilon\delta}\right))=0.$

We pass to the limit as $\varepsilon$ tends to zero in this last equality and we get the desired result.

In order to prove (4.3) let us take $\Psi\in\mathcal{D}\left(\Omega\right)$, $\psi\in H_{per}^{1}(Y)$ and define $\psi_{\varepsilon}(x)=\psi\left(\frac{x}{\varepsilon}\right).$

We have

[TABLE]

By applying the unfolding at the scale $\varepsilon$ we get

[TABLE]

We pass to the limit as $\varepsilon$ tends to zero and we get

[TABLE]

An integration by parts in the domain $Y$ gives

[TABLE]

and this last equality implies (4.3).

In order to prove relation (4.4), let us take $\Psi\in\mathcal{D}\left(\Omega\right).$

We have

[TABLE]

By applying the unfolding at scale $\varepsilon$ and then at scale $\delta$ we get

[TABLE]

We pass to the limit as $\varepsilon$ tends to zero and we get

[TABLE]

which implies (4.4).

Relation (4.5) is a consequence of the following assertions:

[TABLE]

together with the linearity and continuity of the normal trace application from the space $H\left(\operatorname{div},\Omega\right)=\left\{\varphi\in\left(L^{2}\left(\Omega\right)\right)^{n}:\operatorname{div}\varphi\in L^{2}\left(\Omega\right)\right\}$ into $H^{-1/2}\left(\partial\Omega\right).$

By choosing particular test functions in relations (4.2) and (4.3) we obtain relations (4.6) and (4.7) respectively.

Relation (4.8) is a consequence of relations (4.2) and (4.6).

Relation (4.9) is a consequence of relations (4.3) and (4.7).

By applying now the unfolding operator to the inequality (3.4), we get

[TABLE]

In order to pass to the limit in relation (4.10), we consider a test function $v=v_{\varepsilon\delta}$ of the form:

[TABLE]

where $\Psi\in\mathcal{D}\left(\Omega\right)$, $\psi\in\mathcal{D}(Y^{\ast}),$ $\phi\in(H_{per,0}^{1}(Z^{\ast}))^{n},$ $\operatorname{div}_{z}\phi=0.$

We have

[TABLE]

By using this test function we get for the first term in relation (4.10):

[TABLE]

According to the general convergence results for the unfolding we have that the first and second terms tend to zero and the third one to the following limit:

[TABLE]

By using now the fact that the function $B(\varphi)=|\varphi|^{2}$ is proper convex continuous, we have for the fourth term

[TABLE]

In order to pass to the limit in the non linear terms, let us first remark that for a function $v$ in $\left(H^{1}(\Omega)\right)^{n}$ we have

[TABLE]

Indeed, according to a result in [4], we know that

[TABLE]

and following the same ideas we can prove that for a function $w\in L^{2}(\Omega,H^{1}(Y))$ we have

[TABLE]

This implies

[TABLE]

In order to pass to the limit in the first non linear term, by using the previous identity for the function $v_{\varepsilon\delta}$ given by (4.11), we have

[TABLE]

Passing to the limit as $\varepsilon\rightarrow 0,$ we have that

[TABLE]

and

[TABLE]

Moreover, $T_{\delta}(T_{\varepsilon}(\Psi))\rightarrow\Psi$ strongly in $L^{2}\left(\Omega\times Y^{\ast}\times Z^{\ast}\right)$ and so

[TABLE]

Then

[TABLE]

In order to pass to the limit in the second non linear term, we use identity (4.13) for the function $u_{\varepsilon\delta}$ and the fact that the function $E(\varphi)=|\varphi|$ is proper convex continuous. We then deduce:

[TABLE]

Moreover,

[TABLE]

We consider now the term ${{\displaystyle\int_{\Omega\times Y^{\ast}\times Z^{\ast}}}T_{\delta}}\left({T_{\varepsilon}\left(\widetilde{p}_{\varepsilon\delta}\right)}\right)T_{\delta}\left({T_{\varepsilon}\left(\operatorname{div}\left(v-u_{\varepsilon\delta}\right)\right)}\right){dxdydz}$. Using $\operatorname{div}_{x}u_{\varepsilon\delta}=0$ we obtain:

[TABLE]

Passing to the limit as $\varepsilon$ tends to zero and then using (4.2) we obtain

[TABLE]

Combining now all the previous convergences we finally get

[TABLE]

for every $\varPsi\in L^{2}(\Omega,\mathcal{V})$ and by localizing we obtain (4.1).

We notice that the function $\widehat{u}$ which verifies (4.1) is the unique solution of the problem

[TABLE]

for all $\widehat{v}\in\mathcal{V}$.

The non unique function $\widehat{p}$ corresponding to the pressure is then recovered by adapting to our case the ideas in [14]. $\blacksquare$

5 Interpretation of the limit problem

The limit problem (4.1) from Theorem 12 can be interpreted as a non linear Darcy law. In order to derive this result we follow the ideas in Lions and Sanchez-Palencia [14] for the study of the Bingham flow in a classical porous medium.

Let $\lambda\in{\mathbb{R}}^{n}$, $v\in\mathcal{V}$ and define

[TABLE]

Denote $\chi(\lambda)=\chi(y,z;\lambda)$ the unique solution of the following variational inequality: Find $\chi(\lambda)\in\mathcal{V}$ such that

[TABLE]

for every $\Phi\in\mathcal{V}$.

Then we deduce from (4.1) and (5.1) that

[TABLE]

Relations (4.2) and (4.5) imply

[TABLE]

and so the pressure $\widehat{p}$ verifies

[TABLE]

Let us now define

[TABLE]

which is a function from ${\mathbb{R}}^{n}$ into ${\mathbb{R}}^{n}$. Then relation (5.2) reads

[TABLE]

Defining the velocity of filtration by

[TABLE]

we obtain the non linear Darcy law

[TABLE]

where in the right-hand side we have the non linear vectorial function $\mathcal{K}$.

Moreover, according to (4.2) and (4.5), function $u^{0}$ verifies

[TABLE]

[TABLE]

Let us notice that according to Theorem 12 we have

[TABLE]

and

[TABLE]

This clearly shows that (5.3) is the problem verified by the limits of the sequences $u_{\varepsilon\delta}$ and $\widetilde{p}_{\varepsilon\delta}$, solutions of (3.4).

For seek of completness, we recall below the result obtained for the homogenization of the Stokes flow in our porous medium and whose limit is a linear Darcy law. This problem was first studied by Lions in [13] with the method of asymptotic expansions. The justification of the convergence result is done by Bunoiu and Saint Jean Paulin in [7], where the three-scale convergence method introduced by G. Allaire and M. Briane in [2] is used.

The Stokes flow can be seen as a particular case of the Bingham flow and it corresponds to the value zero for the parameter $g$ in the constitutive law. Indeed, when $g$ equals zero, relation (3.1) becomes

[TABLE]

This particular case corresponds to a Newtonian fluid, which satisfies the Stokes system:

[TABLE]

In this case, the unique solution $(u_{\varepsilon\delta},p_{\varepsilon\delta})\in V(\Omega_{\varepsilon\delta})\times L_{0}^{2}(\Omega_{\varepsilon\delta})$ of the Stokes problem satisfies

[TABLE]

Convergence results from Theorem 12 for $u_{\varepsilon\delta}$ and $p_{\varepsilon\delta},$ as far as relations (4.2)-(4.9), hold true. The only difference is the limit problem (4.1) which in this case reads in the simpler way

[TABLE]

for every $\Phi\in\mathcal{V}.$

The linearity of this problem now implies

[TABLE]

where the entries $\chi_{ij},$ $i,j=1,...,n$, of the matrix $\chi$ are the solutions of the following local problems defined in the domain $Y^{\ast}\times Z^{\ast}$:

Find $\chi_{i}\in\mathcal{V}$ such that

[TABLE]

for every $w\in\mathcal{V},$ where $e_{i}$ is the $i$-th unit vector of the canonical base in ${\mathbb{R}}^{n}.$

In this case, the permeability tensor $K$ is defined as the the matrix those entries are

[TABLE]

which is linked for every fixed $i=1,...,n$ to the components of the velocity of filtration via the equality

[TABLE]

where we sum over $j$ between $1$ and $n$.

This is the linear Darcy law for our porous medium, which can be also written as

[TABLE]

where in the right-hand side we multiply a $n\times n$ matrix with a vector belonging to ${\mathbb{R}}^{n}$.

We observe that the linear Darcy law can be seeen as a particular case of the non linear one. Indeed, it is obtained when the function $\mathcal{K}$ of $\lambda$ is linear and so $\mathcal{K}(\lambda)=K\lambda$, where $K$ is a $n\times n$ matrix.