# Analysis of Subgrid Stabilization Method for Stokes-Darcy Problems

Authors: Kamel Nafa
Journal: Advances in Numerical Analysis
 (2016)

---

## Abstract

A number of techniques, used as remedy to the instability of the Galerkin finite element formulation for Stokes like problems, are found in the literature. In this work we consider a coupled Stokes-Darcy problem, where in one part of the domain the fluid motion is described by Stokes equations and for the other part the fluid is in a porous medium and described by Darcy law and the conservation of mass. Such systems can be discretized by heterogeneous mixed finite elements in the two parts. A better method, from a computational point of view, consists in using a unified approach on both subdomains. Here, the coupled Stokes-Darcy problem is analyzed using equal-order velocity and pressure approximation combined with subgrid stabilization. We prove that the obtained finite element solution is stable and converges to the classical solution with optimal rates for both velocity and pressure.


---

## Body

## 1. Introduction

The transport of substances between surface water and groundwater has attracted a lot of interest into the coupling of viscous flows and porous media flows [1–5]. In this work we consider coupled problems in fluid dynamics where the fluid in one part of the domain is described by the Stokes equations and in the other, porous media part, by the Darcy equation and mass conservation. Velocity and pressure on these two parts are mutually coupled by interface conditions derived in [6]. Such systems can be discretized by heterogeneous finite elements as analyzed by Layton et al. [1]. In more recent works, unified approaches become more popular. For instance, discontinuous Galerkin methods were analyzed by Girault and Rivière [3], mixed methods by Karper et al. [4], and local pressure gradient stabilized methods by Braack and Nafa [7].In this work, we consider theL
2-formulation of the coupled Stokes-Darcy problem as in [4], but we discretize by equal-order finite elements and use subgrid method and grad-div term to stabilize the pressure and control the natural H
1
(
d
i
v
) velocity norm on the Darcy subdomain.


## 2. Formulations of the Stokes-Darcy Coupled Equations

### 2.1. Model Equations

LetΩ
⊂
R
d, d
=
2 or 3, be a bounded region with Lipschitz boundary ∂
Ω. Ω
S and Ω
D are, respectively, the fluid and porous media subdomains of Ω such that Ω
S
∩
Ω
D
=
∅. The subdomains have a common interface Γ
=
Ω
¯
S
∩
Ω
¯
D. We denote by v
=
(
v
S
,
v
D
) the fluid velocity and by p
=
(
p
S
,
p
D
) the fluid pressure, where v
i
=
v
|
Ω
i, p
i
=
p
|
Ω
i, i
=
S
,
D
. The flow in the domain Ω
S is assumed to be of Stokes type and governed by the equations (1)
-
2
ν
div
⁡
D
v
S
+
∇
p
S
=
f
,
i
n
Ω
S
div
⁡
v
S
=
0
,
i
n
Ω
Swith symmetric strain tensor D
(
v
S
)
=
1
/
2
(
∇
v
S
+
∇
v
S
T
), external force f
, and constant viscosity ν
>
0
. In the porous region Ω
D the filtration of an incompressible flow through porous media is described by Darcy equations (2)
K
-
1
v
D
+
∇
p
D
=
f
,
i
n
Ω
D
div
⁡
v
D
=
g
,
i
n
Ω
D
,where the permeability K
=
K
(
x
) is a positive definite symmetric tensor and g denotes an external Darcy force.


### 2.2. Boundary Conditions

OnΓ
S
=
∂
Ω
S
∖
Γ, we prescribe homogeneous Dirichlet conditions for the velocity v
S. (3)
v
S
=
0
,
o
n
Γ
S
.The boundary of Ω
D is split into three parts ∂
Ω
D
=
Γ
∪
Γ
D
,
1
∪
Γ
D
,
2
. We prescribe zero flux on Γ
D
,
1 and a homogeneous Dirichlet condition for the pressure on Γ
D
,
2
.
(4)
v
D
·
n
D
=
0
,
o
n
Γ
D
,
1
p
D
=
0
,
o
n
Γ
D
,
2
,where n
D denotes the outer normal vector on the boundary pointing from Ω
D into Ω
S. This boundary condition ensures a zero mass flux.


### 2.3. The Beavers-Joseph-Saffman Condition

The flows inΩ
S and Ω
D are coupled across the interface Γ. Conditions describing the interaction of the flows are as follows [6, 8]:(i)
The continuity of the normal velocity:(5)
v
S
·
n
S
=
-
v
D
·
n
D
,
o
n
Γ
(ii)
The balance of normal forces:(6)
-
-
p
S
I
+
2
ν
D
v
S
n
S
·
n
S
=
p
D
,
o
n
Γ
(iii)
The Beavers-Joseph-Saffman condition written in terms of the strain tensor:(7)
v
S
·
τ
=
-
2
k
~
α
D
v
S
·
n
S
·
τ
,
wherek
~
=
ν
K
τ
·
τ and α is a dimensionless parameter to be determined experimentally, this condition relating the tangential slip velocity v
S
·
τ to the normal derivative of the tangential velocity component in the Stokes region


## 2.1. Model Equations

LetΩ
⊂
R
d, d
=
2 or 3, be a bounded region with Lipschitz boundary ∂
Ω. Ω
S and Ω
D are, respectively, the fluid and porous media subdomains of Ω such that Ω
S
∩
Ω
D
=
∅. The subdomains have a common interface Γ
=
Ω
¯
S
∩
Ω
¯
D. We denote by v
=
(
v
S
,
v
D
) the fluid velocity and by p
=
(
p
S
,
p
D
) the fluid pressure, where v
i
=
v
|
Ω
i, p
i
=
p
|
Ω
i, i
=
S
,
D
. The flow in the domain Ω
S is assumed to be of Stokes type and governed by the equations (1)
-
2
ν
div
⁡
D
v
S
+
∇
p
S
=
f
,
i
n
Ω
S
div
⁡
v
S
=
0
,
i
n
Ω
Swith symmetric strain tensor D
(
v
S
)
=
1
/
2
(
∇
v
S
+
∇
v
S
T
), external force f
, and constant viscosity ν
>
0
. In the porous region Ω
D the filtration of an incompressible flow through porous media is described by Darcy equations (2)
K
-
1
v
D
+
∇
p
D
=
f
,
i
n
Ω
D
div
⁡
v
D
=
g
,
i
n
Ω
D
,where the permeability K
=
K
(
x
) is a positive definite symmetric tensor and g denotes an external Darcy force.


## 2.2. Boundary Conditions

OnΓ
S
=
∂
Ω
S
∖
Γ, we prescribe homogeneous Dirichlet conditions for the velocity v
S. (3)
v
S
=
0
,
o
n
Γ
S
.The boundary of Ω
D is split into three parts ∂
Ω
D
=
Γ
∪
Γ
D
,
1
∪
Γ
D
,
2
. We prescribe zero flux on Γ
D
,
1 and a homogeneous Dirichlet condition for the pressure on Γ
D
,
2
.
(4)
v
D
·
n
D
=
0
,
o
n
Γ
D
,
1
p
D
=
0
,
o
n
Γ
D
,
2
,where n
D denotes the outer normal vector on the boundary pointing from Ω
D into Ω
S. This boundary condition ensures a zero mass flux.


## 2.3. The Beavers-Joseph-Saffman Condition

The flows inΩ
S and Ω
D are coupled across the interface Γ. Conditions describing the interaction of the flows are as follows [6, 8]:(i)
The continuity of the normal velocity:(5)
v
S
·
n
S
=
-
v
D
·
n
D
,
o
n
Γ
(ii)
The balance of normal forces:(6)
-
-
p
S
I
+
2
ν
D
v
S
n
S
·
n
S
=
p
D
,
o
n
Γ
(iii)
The Beavers-Joseph-Saffman condition written in terms of the strain tensor:(7)
v
S
·
τ
=
-
2
k
~
α
D
v
S
·
n
S
·
τ
,
wherek
~
=
ν
K
τ
·
τ and α is a dimensionless parameter to be determined experimentally, this condition relating the tangential slip velocity v
S
·
τ to the normal derivative of the tangential velocity component in the Stokes region


## 3. Variational Formulation

As variational formulation we consider the so-calledL
2- formulation used by Karper et al. [4] and recently by [9, 10]. We denote (8)
v
,
w
Ω
=
∫
Ω
v
w
d
x
,
v
,
w
∈
L
2
Ω
d
,
v
,
w
Γ
=
∫
Γ
v
w
d
s
,
v
,
w
∈
L
2
Γ
,where L
2
(
Ω
) and H
1
(
Ω
) denote the usual Sobolev spaces.Next, we define the spaces(9)
H
Γ
S
1
Ω
S
=
w
∈
H
1
Ω
S
d
∣
w
=
0
o
n
Γ
S
H
1
d
i
v
,
Ω
D
=
w
∈
L
2
Ω
D
d
∣
div
⁡
w
∈
L
2
Ω
D
,
H
Γ
D
,
1
1
Ω
D
=
w
∈
H
1
d
i
v
,
Ω
D
∣
w
·
n
D
=
0
o
n
Γ
D
.Then, multiplying the Stokes equations (1) by the test functions w
S
∈
H
Γ
S
1
(
Ω
S
), q
S
∈
L
2
(
Ω
S
), respectively, and integrating by part on the domain Ω
S, we obtain (10)
2
ν
D
v
S
,
D
w
S
Ω
S
-
2
ν
D
v
S
n
S
,
w
S
-
p
S
,
div
⁡
w
S
Ω
S
+
p
S
,
w
S
·
n
S
Γ
=
f
,
w
S
Ω
S
,
div
⁡
v
S
,
q
S
Ω
S
=
0
.Using the decomposition w
S
=
w
S
·
n
S
n
S
+
w
S
·
τ
τ, the fluid normal stress condition (6), and the BJS interface condition (7) in (10), we obtain the weak formulation of the Stokes equations: find v
S
∈
H
Γ
S
1
(
Ω
S
), p
S
∈
L
2
(
Ω
S
) such that (11)
2
ν
D
v
S
,
D
w
S
Ω
S
+
ν
α
k
~
v
S
·
τ
,
w
S
·
τ
Γ
-
p
S
,
div
⁡
w
S
Ω
S
+
p
D
,
w
S
·
n
S
Γ
=
f
,
w
S
Ω
S
,
div
⁡
v
S
,
q
S
Ω
S
=
0
,
∀
w
S
∈
H
Γ
S
1
(
Ω
S
), q
S
∈
L
2
(
Ω
S
)
.Similarly, takingδ
>
0 and testing the Darcy equations (2) by w
D
∈
H
Γ
D
,
1
1
(
Ω
D
), q
D
∈
L
2
(
Ω
D
), respectively, together with the weighted grad-div term we obtain the weak formulation of Darcy equations: find v
D
∈
H
Γ
D
1
(
Ω
D
), p
D
∈
H
D
,
2
1
(
Ω
D
) such that (12)
K
-
1
v
D
,
w
D
Ω
D
+
∇
p
D
,
w
D
Ω
D
+
δ
div
⁡
v
D
,
div
⁡
w
D
Ω
D
=
δ
g
,
div
⁡
w
D
Ω
D
,
-
v
D
,
∇
q
D
Ω
D
+
v
D
·
n
D
,
q
D
Γ
=
g
,
q
D
Ω
D
.Summing up (11) and (12) the weak form of the coupled problem is given by the following: find v
S
∈
H
Γ
S
1
(
Ω
S
), p
S
∈
L
2
(
Ω
S
), v
D
∈
H
Γ
D
1
(
Ω
D
), and p
D
∈
L
2
(
Ω
D
) such that(13)
2
ν
D
v
S
,
D
w
S
Ω
S
-
p
S
,
div
⁡
w
S
Ω
S
+
K
-
1
v
D
,
w
D
Ω
D
+
∇
p
D
,
w
D
Ω
D
+
δ
div
⁡
v
D
,
div
⁡
w
D
Ω
D
+
ν
α
k
~
v
S
·
τ
,
w
S
·
τ
Γ
+
p
D
,
w
S
·
n
S
Γ
=
f
,
w
S
Ω
S
+
δ
g
,
div
⁡
w
D
Ω
D
,
div
⁡
v
S
,
q
S
Ω
S
-
v
D
,
∇
q
D
Ω
D
-
v
S
·
n
S
,
q
D
Γ
=
g
,
q
D
Ω
D
.To analyze the weak formulation of the coupled problem we introduce the following spaces(14)
V
=
v
∈
H
d
i
v
,
Ω
∣
v
S
∈
H
1
Ω
S
d
,
v
S
=
0
o
n
Γ
S
,
v
·
n
D
=
0
o
n
Γ
D
,
1
,
Q
=
q
∈
L
2
Ω
∣
p
D
∈
H
1
Ω
D
,
p
=
0
∈
Γ
D
,
2
,
X
=
V
×
Q
.The velocity and pressure spaces V and Q are equipped with the natural norms(15)
v
V
=
∇
v
Ω
S
2
+
v
Ω
D
2
+
div
⁡
v
Ω
D
2
1
/
2
,
p
Q
=
p
Ω
S
2
+
∇
p
Ω
D
2
1
/
2
.Further, due to the positive definiteness ofK with respect to the L
2
(
Ω
D
) norm ·
Ω
D, there exist positive real numbers k
1 and k
2 such that (16)
k
1
v
Ω
D
2
≤
K
-
1
v
,
v
Ω
D
≤
k
2
v
Ω
D
2
,
∀
v
∈
V
.Next, we define the bilinear forms forv
=
(
v
S
,
v
D
), w
=
(
w
S
,
w
D
) in V and p
=
(
p
S
,
p
D
), q
=
(
q
S
,
q
D
) in Q on the two parts of the domain by (17)
A
S
v
,
p
;
w
,
q
=
2
ν
D
v
S
,
D
w
S
Ω
S
+
ν
α
k
~
v
S
·
τ
,
w
S
·
τ
Γ
-
p
S
,
div
⁡
w
S
Ω
S
+
div
⁡
v
S
,
q
S
Ω
S
,
A
D
v
,
p
;
w
,
q
=
K
-
1
v
D
,
w
D
Ω
D
+
δ
div
⁡
v
D
,
div
⁡
w
D
Ω
D
+
∇
p
D
,
w
D
Ω
D
-
v
D
,
∇
q
D
Ω
D
.Hence, the bilinear form for the coupled problem is the sum of A
S
(
v
,
p
;
w
,
q
), A
D
(
v
,
p
;
w
,
q
), and terms to enforce the continuity of the normal part of the velocities across the interface. (18)
A
v
,
p
;
w
,
q
=
A
S
v
,
p
;
w
,
q
+
A
D
v
,
p
;
w
,
q
+
p
D
,
w
S
·
n
S
Γ
-
q
D
,
v
S
·
n
S
Γ
.Assuming, for simplicity, that f and g are extended by zero to the whole domain, the variational formulation of the coupled Stokes-Darcy system in compact form reads as follows: find (
v
,
p
)
∈
V
×
Q solution of(19)
A
v
,
p
;
w
,
q
=
F
w
,
q
,
∀
w
,
q
∈
V
×
Q
,with(20)
F
w
,
q
=
f
,
w
S
Ω
+
g
,
q
D
Ω
+
δ
g
,
div
⁡
w
D
Ω
.It can easily be shown that a sufficiently regular solution (
v
,
p
)
∈
V
×
Q of (19) such that v
S
∈
H
2
(
Ω
S
)
d, v
D
∈
H
1
(
Ω
D
)
d, p
∈
H
1
(
Ω
S
∪
Ω
D
) is also a classical solution of (1) and (2). We note that there is an alternative variational formulation to the one given here called H
(
d
i
v
)-formulation. The latter uses the term -
(
p
,
div
⁡
w
)
Ω
D
+
(
div
⁡
v
,
q
)
Ω
D instead of (
w
,
∇
p
)
Ω
D
-
(
v
,
∇
q
)
Ω
D [4].The existence and uniqueness of the solution of problem (19) follows from Brezzi’s conditions for saddle point problems [11]; namely, (21)
A
v
,
p
;
v
,
p
≥
γ
~
⦀
v
⦀
V
2
,
∀
v
∈
V
,
γ
~
>
0
,
(22)
inf
q
∈
L
2
Ω
S
⁡
sup
v
∈
H
1
Ω
S
d
⁡
div
⁡
v
,
q
Ω
S
∇
v
Ω
S
q
Ω
S
≥
β
S
,
(23)
inf
q
∈
H
1
Ω
D
⁡
sup
v
∈
L
2
Ω
D
d
⁡
-
v
,
∇
q
Ω
D
v
Ω
D
∇
q
Ω
D
≥
β
D
.with positive constants β
S and β
D [7].The following lemma is needed in the analysis below and is a consequence of the continuous inf-sup conditions (23) [10].Lemma 1.
For every(
v
,
p
)
∈
X there is w
∈
V such that w
S
·
n
S
=
0 on Γ, satisfying (24)
A
v
,
p
;
w
,
0
≥
c
2
p
Q
2
-
c
1
v
V
2
,
w
V
≤
c
3
p
Q
,with positive constants c
1, c
2, and c
3.Proof.
Let(
v
,
p
)
∈
X. Then, due to Stokes inf-sup condition there exists w
S
∈
H
1
(
Ω
S
)
d with w
S
=
0 on Γ
S and w
S
·
n
=
0 on Γ such that (25)
-
div
⁡
w
S
,
p
Ω
S
=
p
Ω
S
2
,
∇
w
S
Ω
S
≤
c
S
p
Ω
S
.For the Darcy equation, due to the condition p
=
0 on Γ
D
,
2, there exists w
D
∈
H
1
(
Ω
D
)
d with w
D
·
n
=
0 on Γ
D
,
2 and Γ, such that (26)
-
div
⁡
w
D
,
p
Ω
D
=
∇
p
Ω
D
2
,
∇
w
D
Ω
D
≤
c
D
∇
p
Ω
D
.Define (27)
w
=
w
S
i
n
Ω
S
w
D
i
n
Ω
D
,and then (28)
A
v
,
p
;
w
,
0
=
2
ν
D
v
,
D
w
Ω
S
-
p
,
div
⁡
w
Ω
S
+
K
-
1
v
,
w
Ω
D
+
∇
p
,
w
Ω
D
+
δ
div
⁡
v
,
div
⁡
w
Ω
D
≥
-
2
ν
D
v
Ω
S
D
w
Ω
S
+
p
Ω
S
2
-
k
2
v
Ω
D
w
D
Ω
D
+
∇
p
Ω
D
2
-
δ
div
⁡
v
Ω
D
div
⁡
w
Ω
D
≥
-
2
ν
∇
v
Ω
S
∇
w
Ω
S
+
p
Ω
S
2
-
k
2
v
Ω
D
w
D
Ω
D
+
∇
p
Ω
D
2
-
δ
div
⁡
v
Ω
D
∇
w
Ω
D
≥
-
2
ν
c
S
∇
v
Ω
S
p
Ω
S
+
p
Ω
S
2
-
c
p
c
D
k
2
v
Ω
D
∇
p
Ω
D
+
∇
p
Ω
D
2
-
δ
c
D
div
⁡
v
Ω
D
∇
p
Ω
D
,where c
p denote the Poincaré constant.
Then, using Young’s inequality we obtain(29)
A
v
,
p
;
w
,
0
≥
-
ν
c
S
ε
1
∇
v
Ω
S
2
+
1
-
ν
c
S
ε
1
p
Ω
S
2
-
c
p
c
D
k
2
2
ε
2
v
Ω
D
2
+
1
-
c
p
c
D
k
2
ε
2
2
-
δ
c
D
ε
3
2
∇
p
Ω
D
2
-
δ
c
D
2
ε
3
div
⁡
v
Ω
D
2
.
Choosingε
1, ε
2, ε
3 positive constants such that (30)
ε
1
<
1
ν
c
S
,
ε
2
<
2
c
p
c
D
k
2
,
ε
3
<
2
-
c
p
c
D
k
2
ε
2
δ
c
D
,we obtain the required result (31)
A
v
,
p
;
w
,
0
≥
c
2
p
Q
2
-
c
1
v
V
2
,where (32)
c
1
=
max
⁡
ν
c
S
ε
1
,
c
p
c
D
k
2
2
ε
2
,
δ
c
D
2
ε
3
,
c
2
=
min
⁡
1
-
ν
c
S
ε
1
,
1
-
c
p
c
D
k
2
ε
2
2
-
δ
c
D
ε
3
2
.In addition, we also have(33)
w
V
2
=
∇
w
Ω
S
2
+
w
Ω
D
2
+
div
⁡
w
Ω
D
2
≤
c
S
2
p
Ω
S
2
+
c
D
2
c
p
2
+
1
∇
p
Ω
D
2
≤
c
3
p
Q
,where c
3
2
=
m
a
x
c
S
2
,
c
D
2
(
c
p
2
+
1
).


## 4. Finite Element Discretization

LetT
h be a shape-regular partition of quadrilaterals for d
=
2 or hexahedra for d
=
3 [12, 13]. The diameter of element T
∈
T
h will be denoted by h
T and the global mesh size is defined by h
≔
m
a
x
{
h
T
:
T
i
n
T
h
}. Let T
^
≔
(
-
1
;
1
)
d be the reference element, F
T the mapping from T
^ to element T, and Q
r
(
T
^
) the space of all polynomials on T
^ with maximal degree r
≥
0 in each coordinate. We assume that the mesh T
h is obtained from a coarser mesh T
2
h by global refinement. Hence, T
2
h consists of patches of elements of T
h
. We define the finite element space (34)
X
h
r
≔
v
∈
C
Ω
S
∪
C
Ω
D
:
v
T
∘
F
T
i
n
Q
r
T
^
,
∀
T
∈
T
h
.For the discrete spaces V
h and Q
h we use the equal-order finite element functions that are continuous in Ω
S and Ω
D and piecewise polynomials of degree r
≥
1. (35)
V
h
=
X
h
r
d
∩
V
,
Q
h
=
X
h
r
∩
Q
∩
H
1
Ω
.We define theScott-Zhang interpolation operator which preserves the boundary condition [13], as j
r
h
:
H
1
(
Ω
)
→
X
h
r with stability and interpolation properties, respectively, as(36)
∇
j
r
h
ϕ
Ω
≤
c
s
ϕ
1
,
Ω
,
ϕ
∈
H
1
Ω
.
(37)
ϕ
-
j
r
h
ϕ
m
,
Ω
≤
c
i
h
r
+
1
-
m
ϕ
r
+
1
,
Ω
,
ϕ
∈
H
r
+
1
Ω
,
m
=
0
or
1
,where c
i, c
s are positive constants.We will also use the inverse inequality(38)
∑
T
∈
T
h
h
T
2
∇
ϕ
T
1
/
2
≤
c
I
ϕ
Ω
,
∀
ϕ
∈
H
1
Ω
.Similarly, for vector functions we define the interpolation operator(39)
j
r
h
:
H
1
Ω
d
⟶
X
h
r
d
,with interpolation and stability properties as above.It is known that the standard Galerkin discretizations of the Darcy system are not stable for equal-order elements. This instability stems from the violation of the discrete analogue on to the inf-sup condition. One possibility to circumvent this condition is to work with a modified bilinear formA
h
(
·
;
·
) by adding a stabilization term S
h
(
·
;
·
); that is,(40)
A
h
v
h
,
p
h
;
w
,
q
=
A
v
h
,
p
h
;
w
,
q
+
S
h
p
h
;
q
,such that the stabilized discrete problem reads (41)
A
h
v
h
,
p
h
;
w
,
q
=
F
w
,
q
∀
w
,
q
∈
V
h
×
Q
h
.Unlike in [10] where a combination of a generalized mini element and local projection (LPS) is analyzed and in [14] where a method based on two local Gauss integrals for the Stokes equations is used, here we will analyze the problem using a subgrid method [12, 15, 16].For this method the filter, with respect to the global Lagrange interpolantI
2
h, onto a coarser mesh T
2
h is used. Defining κ
2
h
=
I
-
I
2
h the subgrid stabilization term reads(42)
S
h
p
h
;
q
=
∑
M
∈
T
2
h
h
M
γ
∇
κ
2
h
p
h
,
∇
κ
2
h
q
M
,
r
≥
1
,where γ is patchwise constant.A more attractive method from the computational point is obtained using only the fine mesh with smaller stencil. Definingκ
h
=
I
-
I
h the subgrid stabilization term reads (43)
S
h
p
h
;
q
=
∑
K
∈
T
h
h
K
γ
∇
κ
h
p
h
,
∇
κ
h
q
K
,
r
≥
2
.Next, we prove the stability of the discrete coupled Stokes-Darcy problem with respect to the norm(44)
⦀
v
,
p
⦀
h
=
v
V
2
+
p
Q
2
+
S
h
p
;
p
1
/
2
.


## 5. Stability

Theorem 2.
LetT
h be a quasi-regular partition [13]. Then, the following discrete inf-sup condition holds for some positive constant β
~ independent of the mesh size h. (45)
inf
v
h
,
p
h
∈
V
h
×
Q
h
∖
0
,
0
⁡
sup
w
h
,
q
h
∈
V
h
×
Q
h
∖
0
,
0
⁡
A
v
h
,
p
h
;
w
h
,
q
h
⦀
v
h
,
p
h
⦀
h
⦀
w
h
,
q
h
⦀
h
≥
β
~
.Proof.
First, let(
v
h
,
p
h
)
∈
V
h
×
Q
h, and then the diagonal testing combined with Korn’s inequality and the positivity of K
-
1 give (46)
A
h
v
h
,
p
h
;
v
h
,
p
h
=
A
v
h
,
p
h
;
v
h
,
p
h
+
S
h
p
h
;
p
h
≥
α
~
v
V
2
+
S
h
p
h
;
p
h
.In addition, let w be as in Lemma 1, corresponding to (
v
h
,
p
h
)
∈
V
h
×
Q
h, and set z
=
j
r
h
w
-
w. Then,(47)
A
v
h
,
p
h
;
j
r
h
w
,
0
=
A
v
h
,
p
h
;
w
,
0
+
A
v
h
,
p
h
;
z
,
0
≥
c
2
p
h
Q
2
-
c
1
v
h
V
2
+
A
S
v
h
,
p
h
;
z
,
0
+
A
D
v
h
,
p
h
;
z
,
0
.Next, we estimate A
S
(
v
h
,
p
h
;
z
,
0
) and A
D
(
v
h
,
p
h
;
z
,
0
) as follows:(48)
A
S
v
h
,
p
h
;
z
,
0
=
2
ν
D
v
h
,
D
z
Ω
S
+
∇
p
h
,
z
Ω
S
+
ν
α
k
~
v
h
S
·
τ
,
z
S
·
τ
Γ
,where the first two terms are bounded using Cauchy inequality together with the interpolation, stability, and inverse inequalities(49)
ν
D
v
h
,
D
z
Ω
S
≤
ν
D
v
h
Ω
S
D
z
Ω
S
≤
ν
v
h
V
∇
z
Ω
S
≤
ν
c
i
v
h
V
∇
w
Ω
S
≤
ν
c
3
c
i
v
h
V
p
h
Q
,
∇
p
h
,
z
Ω
S
≤
∑
T
∈
T
h
,
T
⊂
Ω
S
h
T
-
2
z
T
2
1
/
2
∑
T
∈
T
h
,
T
⊂
Ω
S
h
T
2
∇
p
h
T
2
1
/
2
≤
∑
T
∈
T
h
,
T
⊂
Ω
S
h
T
-
2
h
T
2
r
∇
w
T
2
1
/
2
c
I
p
h
Ω
S
≤
c
c
i
c
I
∇
w
Ω
S
p
h
Ω
S
≤
c
c
i
c
I
c
3
p
h
Q
2
.The boundary term is bounded using the trace theorem and the H
1- stability by(50)
ν
α
k
~
v
h
S
·
τ
,
z
S
·
τ
Γ
≤
c
Γ
2
ν
α
k
~
v
h
V
∇
z
Ω
S
≤
c
Γ
2
c
s
c
3
ν
α
k
~
v
h
V
p
h
Q
.Hence, by Young inequality with (51)
ϵ
1
=
c
2
8
ν
c
i
c
3
,
ϵ
2
=
c
2
k
~
8
ν
α
c
Γ
2
c
s
c
3we obtain (52)
A
S
v
h
,
p
h
;
z
,
0
≤
c
2
8
c
4
p
h
Q
2
+
c
4
v
h
V
2
,where c
4
=
(
4
(
ν
c
3
c
i
)
2
+
0.25
(
c
Γ
2
c
s
c
3
)
2
)
/
c
2.
For the Darcy bilinear form we have(53)
A
D
v
h
,
p
h
;
z
,
0
=
K
-
1
v
h
,
z
Ω
D
+
δ
div
⁡
v
h
,
div
⁡
z
Ω
D
+
∇
p
h
,
z
Ω
D
=
K
-
1
v
h
,
z
Ω
D
+
δ
div
⁡
v
h
,
div
⁡
z
Ω
D
+
∇
p
h
-
κ
2
h
p
h
,
z
Ω
D
+
∇
κ
2
h
p
h
,
z
Ω
D
≤
K
-
1
v
h
Ω
D
z
Ω
D
+
δ
div
⁡
v
h
Ω
D
z
Ω
D
+
∇
p
h
-
κ
2
h
p
h
Ω
D
z
Ω
D
+
∇
κ
2
h
p
h
Ω
D
z
Ω
D
≤
k
2
v
h
Ω
D
c
i
w
Ω
D
+
δ
div
⁡
v
h
Ω
D
1
+
c
s
w
Ω
D
+
∇
p
h
-
κ
2
h
p
h
Ω
D
c
i
w
Ω
D
+
c
s
∇
p
h
Ω
D
≤
k
2
c
i
c
3
p
h
Q
+
δ
c
3
1
+
c
s
div
⁡
v
h
Ω
D
p
h
Q
+
c
i
c
3
∇
p
h
-
κ
2
h
p
h
Ω
D
+
c
s
p
h
Q
.Then, by Young inequality and (52) we obtain (54)
A
h
v
h
,
p
h
;
j
r
h
w
,
0
≥
5
c
2
8
p
h
Q
2
-
C
v
h
V
2
+
S
h
p
h
;
p
h
.Scaling j
r
h
w we obtain (55)
A
h
v
h
,
p
h
;
j
r
h
w
,
0
≥
p
h
Q
2
-
C
1
v
h
V
2
+
S
h
p
h
;
p
h
.Choosing (
w
h
,
q
h
)
=
(
v
h
,
p
h
)
+
1
/
1
+
C
1
(
j
r
h
w
,
0
) we obtain (56)
A
h
v
h
,
p
h
;
w
h
,
q
h
≥
v
h
V
2
+
1
1
+
C
1
p
h
Q
2
-
C
1
1
+
C
1
v
h
V
2
=
1
1
+
C
1
v
h
V
2
+
p
h
Q
2
=
1
1
+
C
1
⦀
v
h
,
p
h
⦀
h
2
,
⦀
w
h
,
q
h
⦀
h
≤
⦀
v
h
,
p
h
⦀
h
+
1
1
+
C
1
⦀
j
r
h
w
,
0
⦀
h
≤
⦀
v
h
,
p
h
⦀
h
+
C
2
∇
j
r
h
w
Ω
≤
C
3
⦀
v
h
,
p
h
⦀
hwhich implies the required result (57)
inf
v
h
,
p
h
∈
V
h
×
Q
h
∖
0
⁡
sup
w
h
,
q
h
∈
V
h
×
Q
h
∖
0
⁡
A
h
v
h
,
p
h
;
w
h
,
q
h
⦀
v
h
,
p
h
⦀
h
⦀
w
h
,
q
h
⦀
h
≥
β
~
,
withβ
~
=
C
3
-
1
/
(
1
+
C
1
).


## 6. Error Analysis

Theorem 3.
Assume that the solution(
v
,
p
) of the Stokes-Darcy problem (19) is such that (
v
S
,
p
S
)
∈
V
S
∩
H
r
+
1
(
Ω
S
)
d
×
Q
∩
H
l
+
1
(
Ω
S
), (
v
D
,
p
D
)
∈
V
D
∩
H
r
+
1
(
Ω
D
)
d
×
Q
∩
H
l
+
1
(
Ω
D
), and (
v
h
,
p
h
) is the solution of the stabilized problem (41). Then, the following error estimate holds with constants c
1
,
c
2
,
…
,
c
7 independent of h: (58)
⦀
v
-
v
h
,
p
-
p
h
⦀
h
≤
c
1
ν
+
c
2
2
h
2
r
v
r
+
1
,
Ω
S
2
+
c
3
h
+
c
4
δ
2
h
2
r
v
r
+
1
,
Ω
D
2
+
c
5
+
c
6
γ
1
/
2
h
1
/
2
+
c
7
h
2
h
2
l
p
l
+
1
,
Ω
S
2
+
c
5
+
c
6
γ
1
/
2
h
1
/
2
+
c
7
h
2
h
2
l
p
l
+
1
,
Ω
D
2
1
/
2
.Proof.
Using the stability estimate of Theorem3, there exists (
w
h
,
q
h
)
∈
V
h
×
Q
h, with ⦀
(
w
h
,
q
h
)
⦀
h
≤
C
~ satisfying (59)
⦀
j
r
h
v
-
v
h
,
j
l
h
p
-
p
h
⦀
h
≤
1
β
~
A
h
j
r
h
v
-
v
h
,
j
l
h
p
-
p
;
w
h
,
q
h
⦀
w
h
,
q
h
⦀
h
≤
1
β
~
A
h
v
-
v
h
,
p
-
p
h
;
w
h
,
q
h
⦀
w
h
,
q
h
⦀
h
+
1
β
~
A
h
j
r
h
v
-
v
,
j
l
h
p
-
p
;
w
h
,
q
h
⦀
w
h
,
q
h
⦀
h
.Then, by Galerkin orthogonality property, the first term of (59) is bounded by (60)
A
h
v
-
v
h
,
p
-
p
h
;
w
h
,
q
h
⦀
w
h
,
q
h
⦀
h
=
S
h
p
;
q
h
⦀
w
h
,
q
h
⦀
h
≤
S
h
p
;
p
1
/
2
S
h
q
h
;
q
h
1
/
2
⦀
w
h
,
q
h
⦀
h
≤
S
h
p
;
p
1
/
2
.Hence, the approximation properties of κ
2
h and κ
h imply (61)
1
β
~
A
h
v
-
v
h
,
p
-
p
h
;
w
h
,
q
h
⦀
w
h
,
q
h
⦀
h
≤
1
β
~
γ
∇
κ
2
h
p
Ω
∇
κ
2
h
p
Ω
≤
c
1
β
~
-
1
γ
1
/
2
h
l
+
1
/
2
p
l
+
1
,
Ω
.To estimate the second term of (59) we consider separately each individual term of the bilinear form 1
/
β
~
A
h
(
j
r
h
v
-
v
,
j
l
h
p
-
p
;
w
h
,
q
h
).
Next, Cauchy schwarz and Poincaré inequality for the boundary terms imply(62)
1
β
~
A
S
j
r
h
v
-
v
,
j
l
h
p
-
p
;
w
h
,
q
h
≤
β
~
-
1
ν
∇
j
r
h
v
-
v
Ω
S
∇
w
h
Ω
S
+
j
l
h
p
-
p
Ω
S
∇
w
h
Ω
S
∇
j
r
h
v
-
v
Ω
S
q
h
Ω
S
+
ν
α
c
Γ
2
k
~
∇
j
r
h
v
-
v
Ω
S
∇
w
h
Ω
S
≤
β
~
-
1
c
i
C
~
ν
h
r
v
r
+
1
,
Ω
S
+
h
l
+
1
p
l
,
Ω
S
+
h
r
v
r
+
1
,
Ω
S
+
ν
α
c
Γ
2
k
~
h
r
v
r
+
1
,
Ω
S
,
1
β
~
A
D
j
r
h
v
-
v
,
j
l
h
p
-
p
;
w
h
,
q
h
≤
β
~
-
1
k
2
j
r
h
v
-
v
Ω
D
w
h
Ω
D
+
δ
∇
j
r
h
v
-
v
Ω
D
div
⁡
w
h
Ω
D
+
∇
j
l
h
p
-
p
Ω
D
w
h
Ω
D
+
∇
q
h
Ω
D
j
r
h
v
-
v
Ω
D
≤
β
~
-
1
c
i
C
~
k
2
h
r
+
1
v
r
+
1
,
Ω
D
+
δ
h
r
v
r
+
1
,
Ω
D
+
h
l
p
l
+
1
,
Ω
D
+
h
r
+
1
v
r
+
1
,
Ω
D
.Thus, (63)
⦀
j
r
h
v
-
v
h
,
j
l
h
p
-
p
h
⦀
h
≤
c
~
1
ν
+
c
~
2
h
r
v
r
+
1
,
Ω
S
+
c
~
3
h
+
c
~
4
δ
h
r
v
r
+
1
,
Ω
D
+
c
~
5
+
c
~
6
γ
1
/
2
h
1
/
2
+
c
~
7
h
h
l
p
l
+
1
,
Ω
S
+
c
~
5
+
c
~
6
γ
1
/
2
h
1
/
2
+
c
~
7
h
h
l
p
l
+
1
,
Ω
D
.Squaring the norm and applying Young inequality we obtain (64)
⦀
j
r
h
v
-
v
h
,
j
l
h
p
-
p
h
⦀
h
2
≤
4
c
~
1
ν
+
c
~
2
2
h
2
r
v
r
+
1
,
Ω
S
2
+
4
c
~
3
h
+
c
~
4
δ
2
h
2
r
v
r
+
1
,
Ω
D
2
+
4
c
~
5
+
c
~
6
γ
1
/
2
h
1
/
2
+
c
~
7
h
2
h
2
l
p
l
+
1
,
Ω
S
2
+
4
c
~
5
+
c
~
6
γ
1
/
2
h
1
/
2
+
c
~
7
h
2
h
2
l
p
l
+
1
,
Ω
D
2
.Next, we estimate the interpolation error by (65)
⦀
v
-
j
r
h
v
,
p
-
j
l
h
p
⦀
h
2
=
∇
v
-
j
r
h
v
Ω
S
2
+
v
-
j
r
h
v
Ω
D
2
+
div
⁡
v
-
j
r
h
v
Ω
D
2
+
p
-
j
l
h
p
Ω
S
2
+
∇
p
-
j
l
h
p
Ω
D
2
+
S
h
κ
2
h
p
,
κ
2
h
p
≤
c
i
2
h
2
r
v
r
+
1
,
Ω
S
2
+
c
i
2
h
2
r
h
2
+
1
h
2
r
v
r
+
1
,
Ω
D
2
+
c
~
i
2
h
2
+
γ
h
h
2
l
p
l
+
1
,
Ω
S
2
+
c
~
i
2
+
γ
h
h
2
l
p
l
+
1
,
Ω
D
2
.
Adding the interpolation error (64) to the projection error (65) we obtain the required result (66)
⦀
v
-
v
h
,
p
-
p
h
⦀
h
≤
c
1
ν
+
c
2
2
h
2
r
v
r
+
1
,
Ω
S
2
+
c
3
h
+
c
4
δ
2
h
2
r
v
r
+
1
,
Ω
D
2
+
c
5
+
c
6
γ
1
/
2
h
1
/
2
+
c
7
h
2
h
2
l
p
l
+
1
,
Ω
S
2
+
c
5
+
c
6
γ
1
/
2
h
1
/
2
+
c
7
h
2
h
2
l
p
l
+
1
,
Ω
D
2
1
/
2
.Remark 4.
We note that the analysis above holds true for the triangular subgrid interpolationP
r
-
P
r
-
P
r.Remark 5.
Because of the presence of divergence of the velocity and the gradient of the pressure in the discrete norm, the velocity and pressure solutions areO
(
h
r
) and O
(
h
l
), respectively. So, we expect the L
2-asymptotic rates to be O
(
h
r
+
1
) and O
(
h
l
+
1
).


## 7. Numerical Results

As a test model problem we takeΩ
=
(
0,1
)
×
(
0,1
) and split it into Ω
S
=
(
0,1
/
2
)
×
(
0,1
) and Ω
D
=
(
1
/
2,1
)
×
(
0,1
). The interface boundary is Γ
=
{
0.5
,
y
∣
0
<
y
<
1
}. We take ν
=
1, α
=
1, k
~
=
1, and K
=
I and the right hand sides f, g such that the velocity and pressure solution in the two subdomains are given by (67)
u
S
=
y
4
e
x
,
e
y
cos
⁡
2
x
,
x
,
y
∈
Ω
S
u
D
=
y
4
e
x
,
4
y
3
e
x
,
x
,
y
∈
Ω
D
p
=
y
4
e
x
,
x
,
y
∈
Ω
.Note that for this problem forcing terms are needed to balance the equations; notably additional terms are added to the interface conditions in (6) and (7) as follows: (68)
-
-
p
S
I
+
2
ν
D
v
S
n
S
·
n
S
=
p
D
+
g
1
,
o
n
Γ
,
v
S
·
τ
=
-
2
k
~
α
D
v
S
·
n
S
·
τ
o
n
Γ
,where g
1
=
-
2
y
4
e
x, and g
2
=
e
y
cos
⁡
(
2
x
)
+
4
y
3
e
x
-
2
e
y
sin
⁡
(
2
x
).The problem is solved using aQ
1
-
Q
1 velocity-pressure approximation with a two-level subgrid stabilization on a uniform mesh with δ
=
0.4. Rates of convergence for the velocity and pressure errors for h
=
1
/
8
,
1
/
16
,
1
/
32
,
1
/
64, and 1
/
128 are displayed in Tables 1 and 2.Table 1
Rates of convergence for velocity and pressure solution in the Stokes subdomain.
u
-
u
h
0
,
Ω
S
∇
u
-
u
h
0
,
Ω
S
p
-
p
h
0
,
Ω
S
h = 1
8
—
—
—
h = 1
16
1.9303
1.0284
0.8480
h = 1
32
1.9735
1.0208
0.9149
h = 1
64
1.9890
1.0119
0.9511
h = 1
128
1.9951
1.0055
0.9725Table 2
Rates of convergence for velocity and pressure solution in the Darcy subdomain.
u
-
u
h
0
,
Ω
D
d
i
v
u
-
u
h
0
,
Ω
D
p
-
p
h
0
,
Ω
D
h = 1
8
—
—
—
h = 1
16
0.8813
0.8412
1.0416
h = 1
32
0.9534
0.9235
1.0318
h = 1
64
0.9642
0.9514
1.0167
h = 1
128
0.9857
0.9657
1.0085In Table1, we see clearly that the velocity field in the Stokes subdomain is of second-order accuracy with respect to the L
2-norm and first-order accuracy with respect to H
1-seminorm, and the pressure is of first-order accuracy. In addition, In Table 2, we observe that the velocity field and its divergence are of first-order accuracy in the Darcy subdomain, and the pressure is of first-order accuracy with respect to the L
2-norm. So, clearly these results are in agreement with the theoretical results of the previous section and are comparable to the ones found in [2, 5].



---

*Source: 7389102-2016-09-22.xml*