A New Unified Stabilized Mixed Finite Element Method of the Stokes-Darcy Coupled Problem: Isotropic Discretization ()
1. Introduction
There are many serious problems currently facing the world in which the coupling between groundwater and surface water is important. These include questions such as predicting how pollution discharges into streams, lakes, and rivers its way into the water supply. This coupling is also important in technological applications involving filtration. We refer to the nice overview [1] and the references therein for its physical background, modeling, and standard numerical methods. One important issue in the modeling of the coupled Darcy-Stokes flow is the treatment of the interface condition, where the Stokes fluid meets the porous medium. In this paper, we only consider the so-called Beavers-Joseph-Saffman condition, which was experimentally derived by Beavers and Joseph in [2], modified by Saffman in [3], and later mathematically justified in [4] [5] [6] [7].
It is well known that the discretization of the velocity and the pressure, for both Stokes and Darcy problems and the coupled of them, has to be made in a compatible way in order to avoid instabilities. Since, usually, stable elements for the free fluid flow cannot been successfully applied to the porous medium flow, most of the finite element formulations developed for the Stokes-Darcy coupled problem are based on appropriate combinations of stable elements for the Stokes equations with stable elements for the Darcy equations. In [4] [6] [8] - [25], and in the references therein, we can find a large list of contributions devoted to numerically approximate the solution of this interaction problem, including conforming and nonconforming methods.
There are a lot of papers considering different finite element spaces in each flow region (see, for example, [21] [26] [27] and the references therein). In contrast to this, other articles use the same finite element spaces in both regions by, in general, introducing some penalizing terms (ref. for examples [10] [14] [22] and the references therein).
In [22], a conforming unified finite element has been proposed for the modified coupled Stokes-Darcy problem in a plane domain, which has simple and straightforward implementations. The authors apply the classical Mini-element to the whole coupled Stokes-Darcy coupled problem. An a-priori error analysis is performed with some numerical tests confirming the convergence rates.
In this article, we propose a modification of the Darcy problem which allows us to apply a variant nonconforming finite element to the whole coupled Stokes-Darcy problem. We use a variant nonconforming Crouzeix-Raviart finite element method that has so many advantages for the velocities and piecewise constant for the pressures in both the Stokes and Darcy regions, and apply a stabilization term penalizing the jumps over the element edges of the piecewise continuous velocities. We prove that the formulation satisfies the discrete inf-sup condition, obtaining as a result optimal accuracy with respect to solution regularity. Numerical experiments are also presented, which confirm the excellent stability and optimal performance of our method. The difference between our paper and the reference [22] is that our discretization is nonconforming in both the Stokes domain and Darcy domain (in
,
or 3). As a result, additional terms are included in the priori error analysis that measures the non-conformity of the method. One essential difficulty in choosing the unified discretization is that, the Stokes side velocity is in
while the Darcy side velocity is only in
. Thus, we introduce a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element space (larger than the space
used in [14] ). The choice of
[see (34)] is more natural than the one introduced in [14] since the space
approximates only
and not
, while our a priori error analysis is only valid in this larger space.
The rest of the paper is organized as follows. In Section 2 we present the modified coupled Stokes-Darcy problem in
,
or 3, notations and the weak formulation. Section 3 is devoted to the finite element discretization and the error estimation.
In Section 4, we present the results of numerical experiments to verify the predicted rates of convergence. Finally, we offer our conclusion and the further works in Section 5.
2. Preliminaries and Notation
2.1. Model Problem
We consider the model of a flow in a bounded domain
(
or 3), consisting of a porous medium domain
, where the flow is a Darcy flow, and an open region
, where the flow is governed by the Stokes equations. The two regions are separated by an interface
. Let
,
. Each interface and boundary is assumed to be polygonal (
) or polyhedral (
). We denote by
(resp.
) the unit outward normal vector along
(resp.
). Note that on the interface
, we have
. Figure 1 and Figure 2 give a schematic representation of the geometry.
For any function v defined in
, since its restriction to
or to
could play a different mathematical roles (for instance their traces on
), we will set
and
.
![]()
Figure 1. A sketch of the geometry of the problem (case:
).
![]()
Figure 2. A sketch of the geometry of the problem (case
).
In
, we denote by
the fluid velocity and by p the pressure. The motion of the fluid in
is described by the Stokes equations
(1)
while in the porous medium
, by Darcy’s law
(2)
Here,
is the fluid viscosity,
the deformation rate tensor defined by
and
a symmetric and uniformly positive definite tensor representing the rock permeability and satisfying, for some constants
,
is a term related to body forces and
a source or sink term satisfying the compatibility condition
Finally we consider the following interface conditions on
:
(3)
(4)
(5)
Here, Equation (3) represents mass conservation, Equation (4) the balance of normal forces, and Equation (5) the Beavers-Joseph-Saffman conditions. Moreover,
denotes an orthonormal system of tangent vectors on
,
, and
is a parameter determined by experimental evidence.
Equations (1) to (5) consist of the model of the coupled Stokes and Darcy flows problem that we will study below.
2.2. New Weak Formulation
We begin this subsection by introducing some useful notations. If W is a bounded domain of
and m is a non negative integer, the Sobolev space
is defined in the usual way with the usual norm
and semi-norm
. In particular,
and we write
for
. Similarly we denote by
the
or
inner product. For shortness if W is equal to
, we will drop the index
, while for any
,
,
and
, for
. The space
denotes the closure of
in
. Let
be the space of vector valued functions
with components
in
. The norm and the semi-norm on
are given by
(6)
For a connected open subset of the boundary
, we write
for the
inner product (or duality pairing), that is, for scalar valued functions
,
one defines:
(7)
We also define the special vector-valued functions space
(8)
To give the variational formulation of our coupled problem we define the following two spaces for the velocity and the pressure:
equipped with the norm
(9)
and
(10)
Multiplying the first equation of (1) by a test function
and the second one by
, integrating by parts over
the terms involving
and
, yield the variational form of Stokes equations:
(11)
(12)
Using interface conditions (4) and (5) in (11), we obtain:
(13)
(14)
We apply a similar treatment to the Darcy equations by testing the first equation of (2) with a smooth function
and the second on by
, integrating by parts over
the terms involving
, yield the variational form of Darcy equations:
(15)
(16)
Now, incorporating the first boundary interface condition (3) and taking into account that the vector valued functions in
have (weakly) continuous normal components on
(see [28], Theorem 2.5), the mixed variational formulation of the coupled problem (1)-(5) can be stated as follows [10]: Find
that satisfies
(17)
where the bilinear forms
and
are defined on
and
, respectively, as:
By last, the linear forms L and G are defined as:
Theorem 2.1. If
and
, there exists a unique solution
to the problem (17).
Proof. We will establish the existence and uniqueness of a weak solution by using the classical theory of mixed methods so-called Brezzi conditions (see, e.g., [28], Theorem and Corollary 4.1 in Chapter I).
· It is easy to prove that
and
are continuous. It is also clear that F and G are continuous and bounded. In summary, the triangular inequality and Cauchy-Schwarz inequality lead to:
· Now we define the null space of
i.e.
. From the classical Korn’s inequality ( [8], Section 5), it follows that there is a positive constant C such that:
(18)
Let
. Because
, then there exists unique
such that
. Thus, we obtain:
Hence, the coercivity of the bilinear form
,
follows.
· The second Brezzi condition that we need to verify is the inf-sup condition. Let
, then there exists
such that,
(19)
with the estimation
(20)
Thus, from (19) we have
(21)
Using the estimate (20), we obtain
(22)
Hence, the inf-sup condition,
(23)
follows.
Remark 2.1. Note that if g is of mean zero, (17) directly implies that (1), (2) and (3) hold (the differential equations being understood in the distributional sense), while the interface conditions (4) and (5) are imposed in a weak sense. Also, we observe that the mixed variational formulation of the coupled problem (1)-(5) is equivalent to weak formulation (2.4) (and also (2.5) of [11] ), with the particularity that, in our case, for any
, we have that
.
Now we introduce a modification to the Darcy equation, with the purpose in mind of the development of a unified discretization for the coupled problem, that is, the Stokes and Darcy parts be discretized using the same finite element spaces. The modification that we apply to the Darcy equation follows the idea (same argument) given in [22]. Indeed, we observe that taking the second equation of Darcy’s problem (2) we can write, for any
,
(24)
Then, by adding this equation to the first equation of the variational form in (15), we get:
(25)
(26)
From now on, we work with this modified variational form of Darcy equations.
In the same way that before, incorporating the boundary conditions (3) and remambering that, since
, it was (weakly) continuous normal components on
, the variational form of the modified Stokes-Darcy problem can be written as follows: Find
satisfying
(27)
where the bilinear forms
and
are defined on
,
, respectively, as:
and
By last, the linear forms
and G are defined as:
Then, applying the classical theory of mixed methods it follows the well-posedness of the continuous formulation (27).
Theorem 2.2. There exists a unique
solution to modified formulation (27). In addition, there exists a positive constant
, depending on the continuous inf-sup condition constant for
, the coercivity constant for
and the boundedness constants for
and
, such that:
(28)
We end this section with some notation. In 2D, the curl of a scalar function w is given as usual by
while in 3D, the curl of a vector function
is given as usual by
. Finally, let
be the space of polynomials of total degree not larger than k. In order to avoid excessive use of constants, the abbreviations
and
stand for
and
, respectively, with positive constants independent of
or
.
3. A Priori Error Analysis
3.1. Finite Element Discretization
In this subsection, we will use a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element approximation for the velocity and piecewise constant approximation for the pressure.
Let
be a family of triangulations of
with nondegenerate elements (i.e. triangles for
and tetrahedrons for
). For any
, we denote by
the diameter of T and
the diameter of the largest ball inscribed into T and set
(29)
We assume that the family of triangulations is regular, in the sense that there exists
such that
, for all
. We also assume that the triangulation is conform with respect to the partition of
into
and
, namely each
is either in
or in
(see Figures 3-5).
Let
and
be the corresponding induced triangulations of
and
. For any
, we denote by
(resp.
) the set of its edges (
) or faces (
) (resp. vertices) and set
,
. For
we define
Notice that
can be split up in the form
(30)
where
. Note that
is included in
.
![]()
Figure 4. Example of conforming mesh in 2d.
![]()
Figure 5. Example of nonconforming mesh in 2d.
With every edges
, we associate a unit vector
such that
is orthogonal to E and equals to the unit exterior normal vector to
if
. For any
and any piecewise continuous function
, we denote by
its jump across E in the direction of
:
For
, we set (see Figure 6 in 2D for illustration):
(31)
The triplet
with
is finite element [ [29], Page 83]. The local basis functions are defined by:
(32)
where for each
,
is barycentric coordonates of
.
In classical reference element
, the basis fonctions are given in
by:
(33)
Figures 7-9 represent the surface of these functions in space.
The measure of an element or edge/face is denoted by
and
, respectively.
Based on the above notation, we introduce a variant of the nonconforming Crouzeix-Raviart piecewise linear finite element space (larger than the space
used in [14] )
(34)
and piecewise constant function space
where
is the space of the restrictions to T of all polynomials of degree less than or equal to m. The space
is equipped with the norm
while
![]()
Figure 6.
-nonconforming finite element T in 2d.
the norm on
will be specified later on. The choice of
is more natural than the one introduced in [14] since the space
approximates only
and not
, while our a priori error analysis is only valid in this larger space.
Let us introduce the discrete divergence operator
by
(35)
Then, we can introduce two bilinear forms
and
Then the finite element discretization of (27) is to find
such that
(36)
This is the natural discretization of the modified weak formulation (27) except that the penalizing term
is added. This bilinear form
is defined by following the decomposition (37) of
:
(37)
where
and
Here,
is the length (
) or diameter (
) of E. Note that each element of
only contributes with one jump term in
.
Remark 3.1 The Equation (36) have the matrix representation
where
(resp.
) denote the coefficients of
(resp.
) expanded with respect to a basis for
(resp.
).
We are now able to define the norm on
(see [14] ):
In the sequel, we will denote by
and
various constants independent of h. For the sake of convenience, we will define the bilinear form:
From Hölder’s inequality, we derive the boundedness of
and
:
Lemma 3.1. (Continuity of forms) There holds:
(38)
(39)
(40)
(41)
Lemma 3.2. (Coercivity of Ah) There is an
such that:
(42)
Proof. Let
. We have
We introduce the local space
and for
, we define
with the semi-norm:
(43)
Using Young’s inequality and Green formula, we have:
· Estimate
(
or
). We have by Cauchy-Schwarz inequality:
Also, we have:
(44)
Then,
(45)
Hence we deduce
(46)
· Now we estime the term
. By Cauchy-Schwarz, we obtain:
Thus we deduce the estimation:
(47)
Then,
We apply Korn’s discrete inequality [30] and we get:
(48)
Thus
Hence,
(49)
We have,
(50)
(51)
(52)
The estimates (49), (50), (51) and (52), lead to (42). The proof is complete. o
In order to verify the discrete inf-sup condition, we define the space:
(53)
We define also the Crouzeix-Raviart interpolation operator
by:
(54)
(55)
Lemma 3.3. The operator
is bounded: there is a constant
depending on
,
and N such that
(56)
Proof. The proof is similar to ( [14], Lemma 4.5, Page 2695). o
Then, we have the following result
Lemma 3.4. (Inf-Sup condition) There exists a positive constant
depending on
,
and N such that
(57)
Proof. We use Fortin argument [31], i.e. for each
, we find
such that:
Let
. Then from ( [28], Corollary 2.4, Page 24), there exist vectoriel function
satisfying
(58)
, hence
. We take
and we have:
Thus, we obtain
Using the system (58), we have:
(59)
Also,
(60)
From (59) et (60), we deduce:
(61)
The Inf-Sup condition holds and the proof is complete. o
From Lemma 3.2 and Lemma 3.4 we have the following result:
Theorem 3.1. There exists a unique solution
to the problem (27).
3.2. A convergence Analysis
We now present an a priori analysis of the approximation error: The use of nonconforming finite element leads to
, so the approximation error contains some extra consistency error terms. In fact, the abstract error estimates from [14] give the following result:
Lemma 3.5. Let
be the solution of problem (27) and
be the solution of the discrete problem (36). Then we have
(62)
where
and
are the consistency error terms define by:
(63)
(64)
Note that
, thus
.
For estiming the approximation error, we assume that the solution
of problem (27) satisfies the smoothness assumptions:
Assumption 3.1.
1)
,
,
;
2)
,
,
.
We begin with an estimate for the terms
and
.
Lemma 3.6. (Ref. [14] ) There hold:
(65)
(66)
Finally, let us consider the term
. The smoothness assumption of
implies
, thus
. Clearly,
Thus, we have
(67)
where
In order to evaluate the four face integrals, let us introduce two projections operators in the following.
For any
and
, denote by
the constant space of the restrictions to E and
the projection operator from
on to
such that
(68)
The operator
has the property [32]:
(69)
For any
, we let
be the function in
such that
Using inequality (69), we obtain
(70)
Then we have the following lemma:
Lemma 3.7. (Estimation the four face integrals) There holds:
(71)
(72)
(73)
(74)
Proof.
1) Estimate (71): We begin with an estimate for the first term
. For any face
, there exists at least one element
such that
. Then, from condition (68), Höder’s inequality and inequality (70), it follows that
2) Estimate (72):
We have
, hence
.
(75)
Thus,
Furthermore, summing on
faces, we obtain the estimate:
(76)
3) For the terms
and
, we use the same techniques as in the proof of the bounds for
,
, and we obtain:
The proof is complete. o
From Lemma 3.5, Lemma 3.6 and Lemma 3.7, now we derive the following convergence theorem:
Theorem 3.2 Let the solution
of problem (27) satifies the smoothness assumption (Assumption 3.1). Let
be the solution of the discrete problem (36). Then there exists a positive constant C depending on
and
such that:
(77)
4. Numerical Experiments
In this section we present one test case to verify the predicted rates of convergence. The numerical simulations have been performed on the finite element code FreeFem++ [33] [34] in isotropic coupled mesh of Figure 10.
The solutions have been represented by Mathematica software. For simplicity we choose each domain
,
as the unit square,
, and the permeability tensor
is taken to be the identity. The interface
, is the line
, i.e.
like the show Figure 11.
We consider the application
on the open domain
. In
, we define
and we obtain:
(78)
(79)
We choose quadratic pressure
by
(80)
Thus,
![]()
Figure 10. Isotropic mesh
on coupled domain
.
(81)
The exact solution
satisfies the following condition:
(82)
(83)
and the Beavers-Joseph-Saffman interface conditions on
[
]:
(84)
(85)
(86)
Furthermore, we obtain the right-hand terms
define by
(87)
Thus,
in
leads to
and in
,
is given by:
Figure 12 and Figure 13 give an illustration of isotropic and anisotropic meshes.
We study now the convergence of the speed and the pressure in each subdomain on a triangulation
for the non-conforming finite elements of Crouzeix-Raviart. The results obtained are consistent, as shown in Figures 14-17. The order of convergence in norm
approximately equal to 2 for the speed in each domain
,
and of order 1 for the pressure. Then, we represent on the same triangulation, the curves of isovalours of each component of the speed in
,
(see Figures 18-21).
Finally, we study the structure of the stiffness matrix in
on a uniform triangulation. We find that, when the discretization step h becomes more and more infinitely small, the matrix becomes more and more sparse (cf. Table 1), which partly shows the effectiveness of the method digital used.
The exact solutions
and p are represented by Figures 22-24, while the second members
,
,
and
are represented respectively by Figures 25-28.
![]()
Figure 12. Example of isotropic mesh
in 2d.
![]()
Figure 13. Example of anisotropic mesh in 2d.
![]()
Figure 14. Error for the velocity
in
(log/log plot).
![]()
Figure 15. Error for the pressure
in
(log/log plot).
![]()
Figure 16. Error for the velocity
in
(log/log plot).
![]()
Figure 17. Error for the pressure
in
(log/log plot).
![]()
Figure 18. The isovalue of the first velocity component
in
.
![]()
Figure 19. The isovalue of the second velocity component
in
.
![]()
Figure 20. The isovalue of the first velocity component
in
.
![]()
Figure 21. The isovalue of the second velocity component
in
.
![]()
Table 1. Structure of rigidity Matrix on 10 iterations.
5. Summary
· In this contribution, we have investigated a new mixed finite element method to solve the Stokes-Darcy fluid flow model without introducing any Lagrange multiplier. We have proposed a modification of the Darcy problem which allows us to apply a slight variant nonconforming Crouzeix-Raviart element to
the whole coupled Stokes-Darcy problem. The proposed method is probably the cheapest method for Discontinuous Galerkin (DG) approximation of the coupled system, has optimal accuracy with respect to solution regularity, and has simple and straightforward implementations. Numerical experiments have been also presented, which confirm the excellent stability and accuracy of our method.
· Further works: Further it is well known that an internal layer appears at the interface
as the permeability tensor degenerates, in that case anisotropic meshes have to be used in this layer (see for instance [8] [35] ). Hence we intend to extend our results to such anisotropic meshes.
6. Nomenclatures
·
bounded domain
·
: the porous medium domain
·
·
·
·
(resp.
) the unit outward normal vector along
(resp.
)
·
: the fluid velocity
· p: the fluid pressure
· In 2D, the curl of a scalar function w is given as usual by
· In 3D, the curl of a vector function
is given as usual by
namely,
·
: the space of polynomials of total degree not larger than k
·
: triangulation of
·
: the corresponding induced triangulation of
,
· For any
,
is the diameter of T and
is the diameter of the largest ball inscribed into T
·
and
·
: the set of all the edges or faces of the triangulation
·
: the set of all the edges (
) or faces (
) of a element T
·
·
: the set of all the vertices of a element T
·
· For
,
· For
, we associate a unit vector
such that
is orthogonal to E and equals to the unit exterior normal vector to
· For
,
is the jump across E in the direction of
· In order to avoid excessive use of constants, the abbreviations
and
stand for
and
, respectively, with positive constants independent of
or
.
Acknowledgements
The author thanks Professor Emmanuel Creusé (University of Lille 1, France) for having sent us useful documents and for fruitful discussions concerning the numerical tests. We are thankful to the editor and the anonymous reviewers for many valuable suggestions to improve this paper. This article has been posted since April 05, 2019 according to the link https://arxiv.org/abs/1908.01892.
Data Availability
There are no data underlying the findings in this paper to be shared.
Funding
This work has received no funding.