# Decoupling the Stationary Navier-Stokes-Darcy System with the Beavers-Joseph-Saffman Interface Condition

Authors: Yong Cao; Yuchuan Chu; Xiaoming He; Mingzhen Wei
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

This paper proposes a domain decomposition method for the coupled stationary Navier-Stokes
and Darcy equations with the Beavers-Joseph-Saffman interface condition in order to improve the
efficiency of the finite element method. The physical interface conditions are directly utilized
to construct the boundary conditions on the interface and then decouple the Navier-Stokes and
Darcy equations. Newton iteration will be used to deal with the nonlinear systems. Numerical
results are presented to illustrate the features of the proposed method.


---

## Body

## 1. Introduction

The Stokes-Darcy model has been extensively studied in the recent years due to its wide range of applications in many natural world problems and industrial settings, such as the subsurface flow in karst aquifers, oil flow in vuggy porous media, industrial filtrations, and the interaction between surface and subsurface flows [1–8]. Since the problem domain naturally consists of two different physical subdomains, several different numerical methods have been developed to decouple the Stokes and Darcy equations [6, 9–26]. For other works on the numerical methods and analysis of the Stokes-Darcy model, we refer the readers to [27–45].Recently the more physically valid Navier-Stokes-Darcy model has attracted scientists’ attention, and several coupled finite element methods have been studied for it [46–51]. On the other hand, the advantages of the domain decomposition methods (DDMs) in parallel computation and natural preconditioning have motivated the development of different DDMs for solving the Stoke-Darcy model [6, 10–18, 21, 22]. In this paper, we will develop a domain decomposition method for the Navier-Stokes-Darcy model based on Robin boundary conditions constructed from the interface conditions. This physics-based DDM is different from the traditional ones in the sense that they focus on decomposing different physical domains by directly utilizing the given physical interface conditions.The rest of paper is organized as follows. In Section2, we introduce the Navier-Stokes-Darcy model with the Beavers-Joseph-Saffman interface condition. In Section 3, we recall the coupled weak formulation and the corresponding coupled finite element method for the Navier-Stokes-Darcy model. In Section 4, a parallel domain decomposition method and its finite element discretization are proposed to decouple the Navier-Stokes-Darcy system by using the Robin-type boundary conditions constructed from the physical interface conditions. Finally, in Section 5, we present a numerical example to illustrate the features of the proposed method.


## 2. Stationary Navier-Stokes-Darcy Model

In this section we introduce the following coupled Navier-Stokes-Darcy model on a bounded domainΩ=Ωm∪Ωc⊂ℝd,  (d=2,3); see Figure 1. In the porous media region Ωm, the flow is governed by the Darcy system
(1)u→m=-𝕂∇ϕm,∇·u→m=fm.
Here, u→m is the fluid discharge rate in the porous media, 𝕂 is the hydraulic conductivity tensor, fm is a sink/source term, and ϕm is the hydraulic head defined as z+(pm/ρg), where pm denotes the dynamic pressure, z the height, ρ the density, and g the gravitational acceleration. We will consider the following second-order formulation, which eliminates u→m in the Darcy system:
(2)-∇·(𝕂∇ϕm)=fm.Figure 1
A sketch of the porous median domainΩm, fluid domain Ωc, and the interface Γ.In the fluid regionΩc, the fluid flow is assumed to be governed by the Navier-Stokes equations:
(3)u→c·∇u→c-∇·𝕋(u→c,pc)=f→c,(4)∇·u→c=0,
where u→c is the fluid velocity, pc is the kinematic pressure, f→c is the external body force, ν is the kinematic viscosity of the fluid, 𝕋(u→c,pc)=2ν𝔻(u→c)-pc𝕀 is the stress tensor, and 𝔻(u→c)=(∇u→c+∇Tu→c)/2 is the deformation tensor.LetΓ=Ω¯m∩Ω¯c denote the interface between the fluid and porous media regions. On the interface Γ, we impose the following three interface conditions:
(5)u→c·u→c=-u→m·n→m,(6)-u→c·(𝕋(u→c,pc)·n→c)=g(ϕm-z),(7)-τj·(𝕋(u→c,pc)·n→c)=ανdtrace(∏)τj·u→c,
where n→c and n→m denote the unit outer normal to the fluid and the porous media regions at the interface Γ, respectively, τj(j=1,…,d-1) denote mutually orthogonal unit tangential vectors to the interface Γ, and ∏=𝕂ν/g. The third condition (7) is referred to as the Beavers-Joseph-Saffman (BJS) interface condition [52–55].In this paper, for simplification, we assume that the hydraulic headϕm and the fluid velocity u→c satisfy the homogeneous Dirichlet boundary condition except on Γ, that is, ϕm=0 on the boundary ∂Ωm/Γ and u→c=0 on the boundary ∂Ωc/Γ.


## 3. Coupled Weak Formulation and Finite Element Method

In this section we will recall the coupled weak formulation and the corresponding coupled finite element method for the Navier-Stokes-Darcy model with Beavers-Joseph-Saffman condition. Let(·,·)D denote the L2 inner product on the domain D (D=Ωc or Ωm) and let 〈·,·〉 denote the L2 inner product on the interface Γ or the duality pairing between (H001/2(Γ))′ and H001/2(Γ) [5]. Define the spaces
(8)Xc={v→∈[H1(Ωc)]d∣v→=0on∂ΩcΓ},Qc=L2(Ωc),Xm={ψ∈H1(Ωm)∣ψ=0on∂ΩmΓ},
the bilinear forms
(9)am(ϕm,ψ)=(𝕂∇ϕm,∇ψ)Ωm,ac(u→c,v→)=2ν(𝔻(u→c),𝔻(v→))Ωc,bc(v→,q)=-(∇·v→,q)Ωc,
and the projection onto the tangent space on Γ:
(10)Pτu→=∑j=1d-1(u→·τj)τj.With these notations, the weak formulation of the coupled Navier-Stokes-Darcy model with BJS interface condition is given as follows [46–51]: find (u→c,pc,ϕm)∈Xc×Qc×Xm such that
(11)(u→c·∇u→c,v→)Ωc+ac(u→c,v→)+bc(v→,pc)-bc(u→c,q)+am(ϕm,ψ)+〈gϕm,v→·n→c〉-〈u→c·n→c,ψ〉+ανdtrace(∏)〈Pτu→c,Pτv→〉=(fm,ψ)Ωm+(f→c,v→)Ωc+〈gz,v→·n→c〉,∀(v→,q,ψ)∈Xc×Qc×Xm.Assume that we have in hand regular subdivisions ofΩm and Ωc into finite elements with mesh size h. Then one can define finite element spaces Xmh⊂Xm, Xch⊂Xc and Qch⊂Qc. We assume that Xch and Qch satisfy the inf-sup condition [56, 57]
(12)inf0≠qh∈Qchsup0≠v→h∈Xchbc(v→h,qh)∥v→h∥1∥qh∥0>γ,
where γ>0 is a constant independent of h. This condition is needed in order to ensure that the spatial discretizations of the Navier-Stokes equations used here are stable. See [56, 57] for more details of finite element spaces Xmh,  Xch, and Qch that satisfy (12). One example is the Taylor-Hood element pair that we use in the numerical experiments; for that pair, Xch consists of continuous piecewise quadratic polynomials and Qch consists of continuous piecewise linear polynomials.Then a coupled finite element method with Newton iteration for the coupled Navier-Stokes-Darcy model is given as follows [46]: find (u→c,h,pc,h,ϕm,h)∈Xch×Qch×Xmh in the following procedure.(1)
The initial valueu→c,h0 is chosen.(2)
Form=0,1,2,…,M, solve
(13)(u→c,hm+1·∇u→c,hm,v→h)Ωc+(u→c,hm·∇u→c,hm+1,v→h)Ωc+ac(u→c,hm+1,v→h)+bc(v→h,pc,hm+1)-bc(u→c,hm+1,qh)+am(ϕm,hm+1,ψh)+〈gϕm,hm+1,v→h·n→c〉-〈u→c,hm+1·n→c,ψh〉+ανdtrace(∏)〈Pτu→c,hm+1,Pτv→h〉=(u→c,hm·∇u→c,hm,v→h)Ωc+(fm,ψh)Ωm+(f→c,v→h)Ωc+〈gz,v→h·n→c〉,∀(v→h,qh,ψh)∈Xch×Qch×Xmh.(3)
Setu→c,h=u→c,hm+1,  pc,h=p→c,hm+1, and ϕm,h=ϕm,hM+1.


## 4. Physics-Based Domain Decomposition Method

The coupled finite element method may end up with a huge algebraic system, which combines all parts from the Navier-Stokes equations, Darcy equation, and interface conditions together into one sparse matrix. Hence it is often impractical to directly apply this method to large-scale real world applications. In order to develop a more efficient numerical method in this section, we will directly utilize the three physical interface conditions to construct a physics-based parallel domain decomposition method to decouple the Navier-Stokes and Darcy equations.Let us first consider the following Robin condition for the Darcy system: for a given constantγp>0 and a given function ηp defined on Γ,
(14)γp𝕂∇ϕ^m·n→m+gϕ^m=ηp,onΓ.
Then, the corresponding weak formulation for the Darcy part is given by the following: for ηp∈L2(Γ), find ϕ^m∈Xm such that
(15)am(ϕ^m,ψ)+〈gϕ^mγp,ψ〉=(fm,ψ)Ωm+〈ηpγp,ψ〉,∀ψ∈Xm.Second, we can propose the following two Robin-type conditions for the Navier-Stokes equations: for a given constantγf>0 and given functions ηf and η→fτ defined on Γ,
(16)n→c·(𝕋(u→^c,p^c)·n→c)+γfu→^c·n→c=ηf,onΓ,-Pτ(𝕋(u→^c,pc)·n→c)=ανdtrace(∏)Pτu→^c,onΓ.Then, the corresponding weak formulation for the Navier-Stokes equation is given by the following: forηf∈L2(Γ), find (u→^c,p^c)∈Xc×Qc such that
(17)(u→^c·∇u→^c,v→)Ωc+ac(u→^c,v→)+bc(v→,p^c)-bc(u→^c,q)+γf〈u→^c·n→c,v→·u→c〉+ανdtrace(∏)〈Pτu→^c,Pτv→〉=(f→c,v→)Ωc+〈ηf,v→·n→c〉,∀(v→,q)∈Xc×Qc.Our next step is to show that, for appropriate choices ofγf,  γp,  ηf, and ηp, the solutions of the coupled system (11) are equivalent to the solutions of the decoupled equations (15) and (17), and hence we may solve the latter system instead of the former.Lemma 1.
Let(ϕm,u→c,pc) be the solution of the coupled Navier-Stokes-Darcy system (11) and let (ϕ^m,u→^c,p^c) be the solution of the decoupled Navier-Stokes and Darcy equations (15) and (17) with Robin boundary conditions at the interface. Then, (ϕ^m,u→^c,p^c)=(ϕm,u→c,pc) if and only if γf,  γp,  ηf,  η→fτ, and ηp satisfy the following compatibility conditions:
(18)ηp=γpu→^c·n→c+gϕ^m,(19)ηf=γfu→^c·n→c-gϕ^m+gz.Proof.
Adding (15) and (17) together, we obtain the following: given ηp,ηf∈L2(Γ), find (ϕ^m,u^f,p^c)∈Xm×Xc×Qc such that
(20)(u→^c·∇u→^c,v→)Ωc+ac(u→^c,v→)+bc(v→,p^c)-bc(u→^c,q)+am(ϕ^m,ψ)+γf〈u→^c·n→c,v→·n→c〉+〈gϕ^mγp,ψ〉+ανdtrace(∏)〈Pτu→^c,Pτv→〉=(fm,ψ)Ωm+(f→c,v→)Ωc+〈ηf,v→·n→c〉+〈ηpγp,ψ〉,∀(v→,q,ψ)∈Xm×Xc×Qc.
For the necessity of the lemma, we pickψ=0 and v→ such that Pτv→=0 in (11) and (20); then by subtracting (20) from (11), we get
(21)〈ηf-γfv→f·n→c+gϕm-gz,v→·n→c〉=0,iiiiiiiiiiiiiiiiiiiiiiiii∀v→∈XcwithPτv→=0,
which implies (19). The necessity of (18) can be derived in a similar fashion.
As for the sufficiency of the lemma, by substituting the compatibility conditions (18)-(19), we easily see that (ϕ^m,u→^c,p^c) solves the coupled Navier-Stokes-Darcy system (11), which completes the proof.Now we use the decoupled weak formulation constructed above to propose a physics-based parallel domain decomposition method with Newton iteration as follows.(1)
Initial valuesηp0 and ηf0 are guessed. They may be taken to be zero.(2)
Fork=0,1,2,…, independently solve the Darcy and Navier-Stokes equations constructed above. More precisely, ϕmk∈Xm is computed from
(22)am(ϕmk,ψ)+〈gϕmkγp,ψ〉=〈ηpkγp,ψ〉+(fm,ψ)Ωm,(ϕmk,ψ)+〈gϕmkγp,gϕmkγp,〉=〈ηpkγp,ψ〉+∀ψ∈Xm,andu→ck∈Xc and pck∈Qc are computed from the following Newton iteration.
(i)
Initial valueu→ck,0 is chosen for the Newton iteration. For instance, it may be taken to be u→c0,0=0 and u→ck,0=u→ck-1 for k=1,2,….(ii)
Form=0,1,2,…,M, solve
(23)(u→ck,m+1·∇u→ck,m,v→)Ωc+(u→ck,m·∇u→ck,m+1,v→)Ωc+ac(u→ck,m+1,v→)+bc(v→,pck,m+1)-bc(u→ck,m+1,q)+γf〈u→ck,m+1·n→c,v→·n→c〉+ανdtrace(∏)〈Pτu→ck,m+1,Pτv→〉=(u→ck,m·∇u→ck,m,v→)Ωc+〈ηfk,v→·n→c〉+(f→c,v→)Ωc,iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii∀(v→,q,ψ)∈Xc×Qc×Xm.(iii)
Setu→ck=u→ck,M+1 and pck=pck,M+1.(3)
ηpk+1 and ηfk+1 are updated in the following manner:
(24)ηfk+1=γfγpηpk-(1+γfγp)gϕmk+gz,ηpk+1=-ηfk+(γf+γp)u→ck·n→c+gz.Then the corresponding domain decomposition finite element method is proposed as follows.(1)
Initial valuesηp,h0 and ηf,h0 are guessed. They may be taken to be zero.(2)
Fork=0,1,2,…, independently solve the Darcy and Navier-Stokes equations with the Robin boundary conditions on the interface, which are constructed previously. More precisely, ϕm,hk∈Xmh is computed from
(25)am(ϕm,hk,ψh)+〈gϕm,hkγp,ψh〉=〈ηp,hkγp,ψh〉+(fm,ψh)Ωm,∀ψh∈Xmh,andu→c,hk∈Xch and pc,hk∈Qch are computed from the following Newton iteration.
(i)
Initial valueu→c,hk,0 is chosen for the Newton iteration. For instance, it may be taken to be u→c,h0,0=0 and u→c,hk,0=u→c,hk-1 for k=1,2,….(ii)
Form=0,1,2,…,M, solve
(26)(u→c,hk,m+1·∇u→c,hk,m,v→h)Ωc+(u→c,hk,m·∇u→c,hk,m+1,v→h)Ωc+ac(u→c,hk,m+1,v→h)+bc(v→h,pck,m+1)-bc(u→c,hk,m+1,qh)+γf〈u→c,hk,m+1·n→c,v→h·n→c〉+ανdtrace(∏)〈Pτu→c,hk,m+1,Pτv→h〉=(u→c,hk,m·∇u→c,hk,m,v→h)Ωc+〈ηf,hk,v→h·n→c〉+(f→c,v→h)Ωc,∀(v→h,qh,ψh)∈Xch×Qch×Xmh.(iii)
Setu→c,hk=u→c,hk,m+1  and  pc,hk=pc,hk,M+1.(3)
ηp,hk+1 and ηf,hk+1 are updated in the following manner:
(27)ηf,hk+1=γfγpηp,hk-(1+γfγp)gϕm,hk+gz,ηp,hk+1=-ηf,hk+(γf+γp)u→c,hk·n→c+gz.


## 5. Numerical Example

Example 1.
Consider the model problem (2)–(6) with the BJS interface condition (7) on Ω=[0,π]×[-1,1] with Ωm=[0,π]×[0,1] and Ωc=[0,π]×[-1,0]. Choose (ανd/trace(∏))=1, ν=1, g=1, z=0, and 𝕂=K𝕀, where 𝕀 is the identity matrix and K=1. The boundary condition data functions and the source terms are chosen such that the exact solution is given by
(28)ϕm=(ey-e-y)sin(x)et,u→c=[Kπsin(2πy)cos(x)et,(-2K+Kπ2sin2(πy))sin(x)et]T,pc=0.
We divide Ωm and Ωc into rectangles of height h=1/N and width πh, where N is a positive integer, and then subdivide each rectangle into two triangles by drawing a diagonal. The Taylor-Hood element pair is used for the Navier-Stokes equations, and the quadratic finite element is used for the second-order formulation of the Darcy equation.For the coupled finite element method of the steady Navier-Stokes-Darcy model with BJS interface condition, Table1 provides errors for different choices of h. Using linear regression, the errors in Table 1 satisfy
(29)∥u→c,h-u→c∥0≈0.714h3.011,|u→c,h-u→c|1≈3.867h1.987,∥pc,h-pc∥0≈5.123h3.129,∥ϕm,h-ϕm∥0≈0.354h2.998,|ϕm,h-ϕm|1≈1.556h1.995.
These rates of convergence are consistent with the approximation capability of the Taylor-Hood element and quadratic element, which is third order with respect to L2 norm of u→c and ϕm, second order with respect to H1 norm of u→c and ϕm, and second order with respect to L2 norms of pc. In particular, the third-order convergence rate of pc observed above, which is better than the approximation capability of the linear element, is mainly due to the special analytic solution p=0.Table 1
Errors of the finite element method for the steady Navier-Stokes-Darcy model with BJS interface condition.
h
∥
u
→
h
-
u
→
∥
0
∥
u
→
h
-
u
→
∥
1
∥
p
h
-
p
∥
0
∥
ϕ
h
-
ϕ
∥
0
|
ϕ
h
-
ϕ
|
1
1
/
8
1.367
×
1
0
-
3
6.147
×
1
0
-
2
8.002
×
1
0
-
3
6.940
×
1
0
-
4
2.452
×
1
0
-
2
1
/
16
1.687
×
1
0
-
4
1.577
×
1
0
-
2
8.559
×
1
0
-
4
8.687
×
1
0
-
5
6.187
×
1
0
-
3
1
/
32
2.086
×
1
0
-
5
3.978
×
1
0
-
3
9.506
×
1
0
-
5
1.089
×
1
0
-
5
1.553
×
1
0
-
3
1
/
64
2.594
×
1
0
-
6
9.974
×
1
0
-
4
1.121
×
1
0
-
5
1.363
×
1
0
-
6
3.890
×
1
0
-
4
1
/
128
3.235
×
1
0
-
7
2.496
×
1
0
-
4
1.363
×
1
0
-
6
1.705
×
1
0
-
7
9.733
×
1
0
-
5For the parallel DDM withν=1,  K=1,  γf=0.3, and h=1/32, Figures 2 and 3 show the L2 errors of hydraulic head, velocity, pressure, and ηf. We can see that the parallel domain decomposition method is convergent for γf≤γp. Moreover, Figures 4 and 5 show that a smaller γf/γp leads to faster convergence.Convergence for the velocity of the free flow (a) and the hydraulic head of the porous medium flow (b) versus the iteration counterm for the parallel DDM with BJS interface condition.
(a)
(b)Convergence for the pressure of the free flow (a) andηf (b) versus the iteration counter m for the parallel DDM with BJS interface condition.
(a)
(b)Geometric convergence rate of the velocity of the free flow (a) and the hydraulic head of the porous medium flow (b) for the parallel DDM with BJS interface condition.
(a)
(b)Geometric convergence rate of the pressure of the free flow (a) andηf (b) versus the iteration counter m for the parallel DDM with BJS interface condition.
(a)
(b)Then Tables2 and 3 list some L2 errors in velocity, hydraulic head, pressure, and ηf for the parallel domain decomposition method with γf=0.3 and γp=1.2. The data in these two tables indicate the geometric convergence rate γf/γp since all the error ratios are less than (γf/γp)4=(1/4)4=0.0625.Table 2
L
2 errors in velocity and hydraulic head for the parallel DDM with BJS interface condition.
L
2 velocity errors
e
(
i
)
/
e
(
i
-
4
)
L
2 hydraulic head errors
e
(
i
)
/
e
(
i
-
4
)
e
(
0
)
2.342
×
1
0
-
2
6.338
×
1
0
-
1
e
(
4
)  (i=4)
1.225
×
1
0
-
3
0.0523
3.337
×
1
0
-
2
0.0527
e
(
8
)  (i=8)
6.450
×
1
0
-
5
0.0527
1.756
×
1
0
-
3
0.0526
e
(
12
)  (i=12)
3.395
×
1
0
-
6
0.0526
9.246
×
1
0
-
5
0.0527
e
(
16
)  (i=16)
1.787
×
1
0
-
7
0.0526
4.868
×
1
0
-
6
0.0527
e
(
20
)  (i=20)
9.409
×
1
0
-
9
0.0527
2.562
×
1
0
-
7
0.0526Table 3
L
2 errors in pressure and ηf for the parallel DDM with BJS interface condition.
L
2 velocity errors
e
(
i
)
/
e
(
i
-
4
)
L
2 hydraulic head errors
e
(
i
)
/
e
(
i
-
4
)
e
(
0
)
7.268
×
1
0
-
1
5.668
×
1
0
-
2
e
(
4
)  (i=4)
3.826
×
1
0
-
2
0.0526
2.752
×
1
0
-
3
0.0486
e
(
8
)  (i=8)
2.014
×
1
0
-
3
0.0526
1.399
×
1
0
-
4
0.0508
e
(
12
)  (i=12)
1.060
×
1
0
-
4
0.0526
7.233
×
1
0
-
6
0.0517
e
(
16
)  (i=16)
5.579
×
1
0
-
6
0.0526
3.767
×
1
0
-
7
0.0521
e
(
20
)  (i=20)
2.937
×
1
0
-
7
0.0526
1.969
×
1
0
-
8
0.0523Finally, for the preconditioning feature of the domain decomposition method, Table4 shows the number of iterations M is independent of the grid size h. Here, we set γS=0.3,  γD=1.2,  ν=1, and K=1. Let ϕhk,  u→hk, and phk denote the finite element solutions of ϕDk,  u→Sk, and pSk at the kth step of the domain decomposition algorithm. The criterion used to stop the iteration, that is, to determine the value M, is ∥u→hk-u→hk-1∥0+∥ϕhk-ϕhk-1∥0+∥phk-phk-1∥0<ε, where the tolerance ε=10-5.Table 4
The iteration counterM versus the grid size h for both the parallel Robin-Robin domain decomposition method with BJS interface condition.
h
1
/
8
1
/
16
1
/
32
1
/
64
M
19
19
19
19


## 6. Conclusions

In this paper, a parallel physics-based domain decomposition method is proposed for the stationary Navier-Stokes-Darcy model with the BJS interface condition. This method is based on the Robin boundary conditions constructed from the three physical interface conditions. Moreover, it is convergent with geometric convergence rates if the relaxation parameter is selected properly. The number of iteration steps is independent of the grid size due to the natural preconditioning advantage of the domain decomposition methods.



---

*Source: 136483-2013-10-02.xml*