Effective interface conditions for the forced infiltration of a viscous fluid into a porous medium using homogenization

Effective interface conditions for the forced infiltration of a viscous fluid into a porous medium using homogenization

T. Carraro thomas.carraro@iwr.uni-heidelberg.de, Institute for Applied Mathematics, Heidelberg University, 69120 Heidelberg, Germany. TC was supported by the German Research Council (DFG) through project “Modellierung, Simulation und Optimierung der Mikrostruktur mischleitender SOFC-Kathoden” (RA 306/17-2).    C. Goll christian.goll@iwr.uni-heidelberg.de, Institute for Applied Mathematics, Heidelberg University, 69120 Heidelberg, Germany.    A. Marciniak-Czochra anna.marciniak@iwr.uni-heidelberg.de, Bioquant, Heidelberg University, 69120 Heidelberg, Germany. AM-C was supported by ERC Starting Grant ”Biostruct” No. 210680 and Emmy Noether Programme of German Research Council (DFG).    A. Mikelić andro.mikelic@univ-lyon1.frCorresponding author Université de Lyon, Lyon, F-69003, France; Université Lyon 1, Institut Camille Jordan, UMR 5208, Bât. Braconnier, 43, Bd du 11 novembre 1918, 69622 Villeurbanne Cedex, France. The research of A.M. was partially supported by the Programme Inter Carnot Fraunhofer from BMBF (Grant 01SF0804) and ANR. Research visits of A.M. to the Heidelberg University were supported in part by the Romberg professorship at IWR, Heidelberg University, 2011–1013.

It is generally accepted that the effective velocity of a viscous flow over a porous bed satisfies the Beavers-Joseph slip law. To the contrary, in the case of a forced infiltration of a viscous fluid into a porous medium the interface law has been a subject of controversy. In this paper, we prove rigorously that the effective interface conditions are: (i) the continuity of the normal effective velocities; (ii) zero Darcy’s pressure and (iii) a given slip velocity. The effective tangential slip velocity is calculated from the boundary layer and depends only on the pore geometry. In the next order of approximation, we derive a pressure slip law. An independent confirmation of the analytical results using direct numerical simulation of the flow at the microscopic level is given, as well.

Keywords Interface conditions, pore scale simulation, pressure slip law, slip velocity, boundary layers, homogenization

1 Introduction

The purpose of this paper is to derive rigorously the interface conditions governing the infiltration of a viscous fluid into a porous medium.

We start from an incompressible 2D flow of a Newtonian fluid penetrating a porous medium. At the pore scale, the flow is described by the stationary Stokes system, both, in the unconstrained fluid part and in the pore space. The upscaling of the Stokes system in a porous medium yields Darcy’s law as the effective momentum equation, valid at every point of the porous medium. The Stokes system and the Darcy equation are very different PDEs and need to be coupled at the interface of the fluid and the porous medium. The resulting system should be an approximation of the starting first principles with error estimate in the term of the dimensionless pore size , being the ratio of the characteristic pore size and the macroscopic domain length.

There is vast literature on modeling interface conditions between a free flow and a porous medium. Most of the references focus on flows which are tangential to the porous medium. In such situation, the free fluid velocity is much larger than the Darcy velocity in the porous medium. The corresponding interface condition is the slip law by Beavers and Joseph. It was deduced from the experiment in [3], then discussed and simplified into a generally used form in [31] and justified through numerical simulations of pore level Navier-Stokes equations in [32], [20] and [7]. A rigorous justification of the slip law by Beavers and Joseph, starting from the pore level first principles, was provided by Jäger and Mikelić in [17], using a combination of homogenization and boundary layers techniques. The slip law is supplemented by the pressure jump law, what was noticed in [18] and rigorously derived in [25]. A corresponding numerical validation by solving the Stokes equation at the pore scale has been recently presented in [7].

Infiltration into a porous medium corresponds to a different situation, because in this case the free fluid velocity and the Darcy velocity are of the same order. We refer to the article by Levy and Sanchez-Palencia [23]. They classify the physical situation as ”Case B: The pressure gradient on the side of the porous body at the interface is normal to it”. In the ”Case B” the pressure gradient in the porous medium is much larger than in the free fluid. Using an order-of-magnitude analysis, in [23] it was concluded that the effective interface conditions have to satisfy


where are the Darcy velocity and the pressure and is the unconfined fluid velocity. Note that the interface conditions (1.1) were obtained for low Reynolds number flows.

In order to couple the Stokes system in the free fluid domain with the Darcy equation in the porous medium, the conditions (1.1) are not sufficient. One more condition is needed. In [23] an intermediate boundary layer was introduced and existence of an effective slip velocity at the interface was postulated. However, the article [23] did not provide the slip velocity. It was limited to a model of macroscopic isotropy, where the slip is equal to zero. Therefore, zero tangential velocity is assumed.

A rigorous mathematical study of the interface conditions between a free fluid and a porous medium was initiated in [16]. Our analysis repose on the boundary layers constructed there. For reviews of the models and techniques we refer to [11] and [19].

We note that in a number of articles devoted to numerical simulations, the porous part was modeled using the Brinkman-extended Darcy law. We refer to [10], [13], [14], [27], [36] and references therein. In such setting, the authors used general interface conditions introduced by Ochoa-Tapia and Whitaker in [29]. They consists of (i) continuity of the velocity and (ii) complex jump relations for the stresses, containing several parameters to be fitted. We recall that the viscosity in the Brinkman equation is not known and the use of it seems to be justified only in the case of a high porosity (see the discussion in [28]). Furthermore, Larson and Higdon undertook a detailed numerical simulation of two configurations (axial and transverse) of a shear flow over a porous medium in [21]. Their conclusion was that a macroscopic model based on Brinkman’s equation gives “reasonable predictions for the rate of decay of the mean velocity for certain simple geometries, but fails for to predict the correct behavior for media anisotropic in the plane normal to the flow direction”. An approach using the thermodynamically constrained averaging theory was presented in [15]. Darcy-Navier-Stokes coupling yields also an increasing interest from the side of numerical analysis, see [22], [30] and [11] and references therein.

In this paper, we provide a rigorous derivation of the filtration equation and the interface condition explained above from the pore scale level description based on first principles. Our derivation follows the general homogenization and boundary layers approach from [16]. The necessary results on boundary layers and very weak solutions to the Stokes system will be recalled in the proofs of the main results.

In our work we have used a finite element method to obtain a numerical confirmation of the analytical results. The numerical study of the convergence rates of the macroscopic problems and effective interface conditions is a difficult task for the reasons explained in the following. The first difficulty is the numerical solution of the microscopic problem used as reference. The geometry of the porous part has to be resolved with high accuracy. In addition, the microscopic solution in the vicinity of the surface of the porous medium has large gradients, that can only be approximated by a boundary layer as shown in this work. The accuracy needed by the resolution of the interface and porous part requires high performance computing. We have reduced the computational costs by considering for our test cases a problem with periodic geometry and periodic boundary conditions. We could thus reduce our computations to one column of inclusions in the porous part. Nevertheless, even in our simplified example problem all the computations must be performed with high accuracy. The reason is that the homogenization errors, especially in the estimates based on correction terms, are small in comparison with numerical errors even for simulations with millions of degrees of freedom. A further difficulty is that to numerically check the estimates we have to solve several auxiliary problems that are coupled. Therefore the numerical precision of one problem influences the precision of the other ones. Due to the complexity of the microscopic flow and the boundary layers, strategies for local mesh adaptivity to reduce the computations of the norms in the estimates are not effective. We could nevertheless apply a goal oriented adaptive method, based on the dual weighted residual (DWR) method [4], to calculate some constants needed for the estimates, increasing the overall accuracy of our numerical tests.

The paper is organized as follows: In Section 2 we formulate the starting microscopic problem and the resulting effective equations. We provide theorems on error estimates of the model approximation. In Section 3 we give a numerical confirmation of the analytical results based on finite element computations. Sections 45 contain the corresponding proofs.

2 Problem setting and main results

2.1 Definition of the geometry

Let and be positive real numbers. We consider a two dimensional periodic porous medium with a periodic arrangement of the pores. The formal description goes along the following lines: First, we define the geometrical structure inside the unit cell . Let (the solid part) be a closed strictly included subset of , and (the fluid part). Now we make a periodic repetition of all over and set , . Obviously, the resulting set is a closed subset of and in an open set in . We suppose that has a boundary of class , which is locally located on one side of their boundary. Obviously, is connected and is not. Now we notice that is covered with a regular mesh of size , each cell being a cube , with . Each cube is homeomorphic to , by linear homeomorphism , being composed of translation and a homothety of ratio .

We define For sufficiently small we consider the set and define

Obviously, . The domains and represent, respectively, the solid and fluid parts of the porous medium . For simplicity, we suppose .

We set , and . Furthermore, let and .


Figure 1: Sketch of the geometry LABEL:sub@afoto the periodicity cell   LABEL:sub@subfig.unit_cell.

2.2 The microscopic model

Having defined the geometrical structure of the porous medium, we precise the flow problem.

We consider the slow viscous incompressible flow of a single fluid through a porous medium. The flow is caused by the fluid injection at the boundary We suppose the no-slip condition at the boundaries of the pores (i.e. the filtration through a rigid porous medium). Then, the flow is described by the following non-dimensional steady Stokes system in (the fluid part of the porous medium ):


Such flow is possible only under the following compatibility condition


With the assumption on the geometry from section 2.1, condition (2.2) and for , and , problem (2.1a)-(2.1d) admits a unique solution , for all .

Our goal is to study behavior of solutions to system (2.1a)-(2.1d), when . In such limit the equations in remain unchanged and in the Stokes system is upscaled to Darcy’s equation posed in . Our contribution is the derivation of the interface condition, linking these two systems.

2.3 The boundary layers and effective coefficients

Let the permeability tensor be given by




Obviously, these problems always admit unique solutions and is a symmetric positive definite matrix (the dimensionless permeability tensor). In addition


In order to formulate the result we need the viscous boundary layer problem connecting free fluid flow and a porous medium flow:

Figure 2: The boundary layer geometry

On the figure, the interface is , the free fluid slab is and the semi-infinite porous slab . The flow region is then . We consider the following problem: Find , , with square-integrable gradients satisfying


By Lax-Milgram’s lemma, there exists a unique satisfying (2.6a)-(2.6e) and , which is unique up to a constant and satisfying (2.6a). After the results from [16], the system (2.6a)-(2.6e) describes a boundary layer, i.e. and stabilize exponentially towards constants, when : There exists and and such that


The case is of special importance. If we suppose the mirror symmetry of the solid obstacle with respect to , then it is easy to prove that is uneven in with respect to the line , and and are even. Consequently, and the permeability tensor is diagonal. Next we see that is uneven in with respect to the line , and and are even. Using formula (2.9) yields in the case of the mirror symmetry of the solid obstacle with respect to .

2.4 The macroscopic model

Now we introduce the effective problem in . It consist of two problems, which are to be solved sequentially. The first problem is posed in and reads:

Find a pressure field which is the periodic in function satisfying


We note that the pressure field is equal to a particular constant, which is equal to zero.

Problem (2.11a)-(2.11b) admits a unique solution . Next, we study the situation in the unconfined fluid domain :

Find a velocity field and a pressure field such that


The constant is given by (2.7) and requires solving problem (2.6a)-(2.6e).

Again, using the compatibility condition (2.2) we obtain easily that problem (2.12a)-(2.12d) has a unique solution .

2.5 The main result

In this section we formulate the approximated model. We expect that the Stokes system remains valid in .

Since , the filtration velocity has to be of order . Therefore, after [1], [19] and [33], the asymptotic behavior of the velocity and pressure fields in the porous part , in the limit , is given by the two-scale expansion

The boundary layers given by (2.6a)-(2.6e) will be used to link the above approximation on with the solution of the Stokes system. With such strategy, at the main order approximation reads


where is the Heaviside function. We will see that .

Theorem 1.

Let be a neighborhood of . Let us suppose the geometry and data smoothness as above and the compatibility condition (2.2). Let be extended to by formula (4.80) of Lipton and Avellaneda. Then we have


Inspection of the proof of theorem (1) shows that we can obtain slightly better estimates by rearranging the term

We obtain

Theorem 2.

Let be a neighborhood of . Let us suppose the geometry and data smoothness as above and the compatibility condition (2.2). Let be extended to by formula (4.80) of Lipton and Avellaneda. Then we have

Remark 3.

We took as correction to the quantity . In fact the better choice would be to take a function satisfying equations (2.11a)-(2.11b), with value on being , instead of zero. Since the order of approximation does not change, we make the simplest possible choice.

If the effective porous medium pressure is , then the requirement that we can only have an normal stress jump on yields


The relation (2.23) indicates presence of an effective pressure slip at the interface . Since is related to the square root of the permeability, in the dimensional formulation, our result compares to the numerical experiments by Sahraoui and Kaviany in [20] and [32]. They found it being small for parallel flows. We find it small but of order of the corrections in the law by Beavers and Joseph in the case of the transverse flow.

Remark 4.

We obtain an expression for the tangential slip velocity. Since it is zero in the isotropic case, we do not confirm formulas like (86), page 2645, from [29] or like formula (31) for oblique flows from [32], which generalize the law by Beavers and Joseph.

In [29], formula (71), page 2643, expresses the continuity of the averaged velocities. By construction, we have the trace continuity for our approximation. Nevertheless, one usually does not keep the boundary layers in the macroscopic model. If we eliminate the boundary layers and all low order terms, the tangential effective velocity at the interface is

(see (2.12d) and from the porous media side

Therefore we find out that there is an effective tangential velocity jump at the interface.

3 Numerical confirmation of the effective interface conditions

This section is dedicated to the numerical confirmation of the analytical results shown above. We solve the problems needed to numerically compare the microscopic with the macroscopic problem by the finite element method (FEM). For the FEM theory we refere to standard literature, e.g., [8] or [5].

For the discretization of the Stokes system we use the Taylor-Hood element, which is inf-sup stable [6], therefore it does not need stabilization terms. In particular, since the homogenization error in some of the proposed estimates is small in comparison with the discretization error even for meshes with a number of elements in the order of millions, we have used higher order finite elements (polynomial of third degree for the velocity components and of second degree for the pressure) to reduce the discretization error.

The flow properties depend on the geometry of the pores. In particular there is a substantial difference between the case with symmetric inclusions with respect to the axis orthogonal to the interface and the case with asymmetric inclusions. We use therefore two different types of inclusion in the porous part, circles and rotated ellipses, i.e. ellipses with the major principal axis non parallel to the flow. The increased accuracy using higher order finite elements in the numerical solutions was necessary, as shown later, especially for the case with symmetric inclusions. The geometries of the unit cells , see figure 3, for these two cases are as follows:

  1. the solid part of the cell is formed by a circle with radius and center .

  2. consists of an ellipse with center and semi-axes and , which are rotated anti-clockwise by .

In addition, since the considered domains have curved boundaries we use cells of the FEM mesh with curved boundaries (a mapping with polynomial of second degree was used for the geometry) to obtain a better approximation.

(a) Circle
(b) Ellipse
Figure 3: Mesh of the fluid part of the unit cell for the two types of inclusions: circles and ellipses

All computations are done using the toolkit DOpElib ([12]) based upon the C++-library deal.II ([2]).

3.1 Numerical setting

In this subsection we describe the setting for the numerical test. To confirm the estimates of Theorem 1 and 2 we have to solve the microscopic problem (2.1) to get and , the macroscopic problems (2.11) and (2.12) to get and , the cell problem (2.4) to calculate the permeability tensor , the velocity vector and pressure , and the boundary layer (2.6) for the velocity and pressure .

To reduce the discretization errors we consider a test case, described below, for which it is easy to derive the exact form of the macroscopic solution. As we will show below, the analytical solution of the macroscopic problem can be expressed in terms of the solution of the cell and boundary layer problems. The discretization error of the macroscopic problem in this case depends on the discretization error of the cell and boundary layer problems and does not imply therefore an additional discretization error.

We consider the following domains and , where the obstacles are either circles or ellipses as described in the subsection above. In our example we consider the in- and outflow condition


in the microscopic problem (2.1)

The macroscopic solution in this setting is


The macroscopic solution depends on the solution of the cell problem though the permeability , see expressions (3.2a) and (3.2d). Furthermore it depends on the solution of the boundary layer though the constant . The macroscopic problems (2.11) and (2.12) are therefore not numerically solved.

The microscopic problem (2.1) is solved with around 10–15 million degrees of freedom, the cell problem uses around 7 million degrees of freedom. The permeability constant has been precisely calculated using the goal oriented strategy for mesh adaptivity described in [7].

In the boundary layer problem, due to the interface condition (2.6c), the velocity as well as the pressure are discontinuous on the interface . Since with the conform finite elements chosen for the discretization the discontinuity cannot be properly approximated, we have decided to transform the problem so that the solution variables are continuous across . The values of and needed to check the estimates are recovered by post-processing. For the numerical solution, as explained in detail in the appendix of our previous work [7], we use a cut-off domain, which is justified by the exponential decaying of the boundary layer solution. The solution of the boundary layer problem is obtained with a mesh of around 4 million degrees of freedom and the constants and are calculated by the goal oriented strategy for mesh adaptivity described in [7] where we have made sure that the cut-off error is smaller than the discretization error. We note that in the computation of we do not use the formula given in (2.10) but the equivalent one


as this proved to be advantageous numerically.

In table 1 the computed constants and for the two different inclusions are listed. As the permeability tensor has for the given cases the form

i.e. it holds and , we state only and . Additionally, we give an estimation of the discretization error.

circular inclusions oval inclusions
0.0199014353519271 0.0122773324576884
0 0.00268891986291451
0 -0.003336740001686
0.025777570627281 -0.004429782196436
Table 1: Computed constants for the computations in the example.

3.2 Numerical results

In this section we present the numerical confirmation of the convergence rates of the homogenization errors (2.152.18) and (2.192.22).

(a) Velocity estimates
(b) Velocity estimate in
(c) Velocity on .
(d) Pressure estimates multiplied by
Figure 4: Convergence results for circles.

For our test we set , and we use a computation of the boundary layer on a cut-off domain ranging from to . This means that to compute the norms we evaluate the terms involving the boundary layer only for with . Outside of this region we assume the difference between the boundary layer components and their respective asymptotic values to be sufficiently small.

In the case of inclusions symmetric in the sense explained above, e.g. circles, the homogenization errors are much smaller than the numerical error even for large epsilon such as as can be observed in figure 4). The lines with markers represent the results of the computations for , the solid lines are reference values for various convergence rates and are plotted only to compare the respective slopes.

The case of circles is shown in figure 4. For the velocity in the fluid part of the domain the estimate (2.15) can be verified. For the better estimate (2.19), that uses correction terms to improve the estimation, the homogenization error is so small that the curve shows only the numerical error, that in our case is only due to the discretization error since the quadrature error and the tolerance of the solver are smaller. In Figure (b)b we can confirm (2.21) only for values of epsilon not bigger than , for the numerical error dominates the homogenization error. In the estimates for circles on the interface (Figure (c)c) we can observe only the numerical error for the same reason explained above. Notice that the error for circles shown in Figure 4 is much smaller than the error for ellipses shown in Figure 5. In addition, we could verify both estimates for the pressure (2.18) and (2.22) as shown in Figure (d)d. Note that the pressure estimates have been scaled multiplying by .

(a) Velocity estimates in
(b) Velocity estimate in
(c) Velocity components on
(d) Pressure estimates multiplied by
Figure 5: Convergence results for elliptical inclusions.

The case of ellipses is shown in figure 5. As it can be observed, all estimates could be numerically verified, since the discretization error in this case was smaller than the homogenization error. Also in this case the pressure estimates have been scaled multiplied by . Note, that we observe for the velocity in the porous domain a convergence rate of 1.5 instead of the predicted first order convergence, see figure (b)b.

In conclusion, we show in figure 6 and figure 7 pictures of the flow for the case . Since we use periodic boundary conditions in the -direction, constant in- and outflow data as well as a periodic geometry, the computations have been performed on a stripe of one column of inclusions to reduce the computational effort. In figure (a)a and (c)c we see streamline plots of the velocity, figure (b)b and (d)d show the corresponding pressures. Both pressures are nearly constant in the fluid part and show then a linear descent to the outflow boundary, similar to the effective pressure (3.2c) and (3.2d).

Figure 7 shows only the values of the tangential velocity component. In the case of circular inclusions (figure (a)a), the velocity is nearly zero throughout the fluid region and shows some oscillations around the mean value zero on the position of the interface. Note that the effective model prescribes here a no slip condition because it holds . In figure (b)b) we see the corresponding solution for oval inclusions. We notice a linear descent from the inflow boundary (which lies in this picture on the left hand side) to the interface, which leads to the slip condition for the tangential velocity component of the effective flow in this case. Both behaviors are predicted from the effective interface condition for this velocity component, see (2.12d).

(a) Streamlines of
(c) Streamlines of
Figure 6: Visualization of the solution to the microscopic problem for . Subfigures LABEL:sub@subfig.vcand LABEL:sub@subfig.pc show the results for circular inclusions, LABEL:sub@subfig.veand LABEL:sub@subfig.pe for elliptical inclusions.
(a) Circular inclusions.
(b) Elliptical inclusions.
Figure 7: Visualization of for .

4 Proof of Theorem 1 via incremental accuracy correction

In the proofs which follow we will frequently use the space


We will follow the strategy from [16], write a variational equation for the errors in velocity and in pressure and reduce the forcing term in several steps. We will frequently use the notation


where is given by (2.4).

4.1 Incremental accuracy correction, the 1st part

Proposition 5.

Let given by (2.11a)-(2.11b), be the solution for (2.12a)-(2.12d) and defined by (4.2). Let be the solution for (2.1a)-(2.1d). Then for every we have


where on .

Proof of proposition 5 We start with the weak formulation corresponding to (2.1a)-(2.1d):


As a first step we eliminate the boundary conditions. The weak formulation corresponding to system (2.12a)-(2.12d) is


Next the weak formulation corresponding to the correction in the pore space is


We observe that difference between (4.4) and (4.5)-(4.1) is equivalent to