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

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.

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 Section 2, 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 Ω=Ω_𝑚∪Ω_𝑐⊂ℝ^𝑑,  (𝑑=2,3); see Figure 1. In the porous media region Ω_𝑚, the flow is governed by the Darcy system 𝑢_𝑚=−𝕂∇𝜙_𝑚,∇⋅𝑢_𝑚=𝑓_𝑚. (1) Here, 𝑢_𝑚 is the fluid discharge rate in the porous media, 𝕂 is the hydraulic conductivity tensor, 𝑓_𝑚 is a sink/source term, and 𝜙_𝑚 is the hydraulic head defined as 𝑧+(𝑝_𝑚/𝜌𝑔), where 𝑝_𝑚 denotes the dynamic pressure, 𝑧 the height, 𝜌 the density, and 𝑔 the gravitational acceleration. We will consider the following second-order formulation, which eliminates 𝑢_𝑚 in the Darcy system: −∇⋅𝕂∇𝜙_𝑚=𝑓_𝑚. (2)

Figure 1 

A sketch of the porous median domain Ω𝑚, fluid domain Ω𝑐, and the interface Γ.

In the fluid region Ω_𝑐, the fluid flow is assumed to be governed by the Navier-Stokes equations: 𝑢_𝑐⋅∇𝑢_𝑐−∇⋅𝕋𝑢_𝑐,𝑝_𝑐=𝑓_𝑐,∇⋅𝑢_𝑐=0, (3)(4) where 𝑢_𝑐 is the fluid velocity, 𝑝_𝑐 is the kinematic pressure, 𝑓_𝑐 is the external body force, 𝜈 is the kinematic viscosity of the fluid, 𝕋(𝑢_𝑐,𝑝_𝑐)=2𝜈𝔻(𝑢_𝑐)−𝑝_𝑐𝕀 is the stress tensor, and 𝔻(𝑢_𝑐)=(∇𝑢_𝑐+∇^𝑇𝑢_𝑐)/2 is the deformation tensor.

Let Γ=Ω_𝑚∩Ω_𝑐 denote the interface between the fluid and porous media regions. On the interface Γ, we impose the following three interface conditions: 𝑢_𝑐⋅𝑢_𝑐=−𝑢_𝑚⋅𝑛_𝑚,−𝑢_𝑐⋅𝕋𝑢_𝑐,𝑝_𝑐⋅𝑛_𝑐=𝑔𝜙_𝑚−𝑧,−𝝉_𝑗⋅𝕋𝑢_𝑐,𝑝_𝑐⋅𝑛_𝑐=𝛼𝜈√𝐝(∏)𝝉trace^𝑗⋅𝑢_𝑐, (5)(6)(7) where 𝑛_𝑐 and 𝑛_𝑚 denote the unit outer normal to the fluid and the porous media regions at the interface Γ, respectively, 𝝉_𝑗(𝑗=1,…,𝑑−1) denote mutually orthogonal unit tangential vectors to the interface Γ, and ∏=𝕂𝜈/𝑔. 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 𝜙_𝑚 and the fluid velocity 𝑢_𝑐 satisfy the homogeneous Dirichlet boundary condition except on Γ, that is, 𝜙_𝑚=0 on the boundary 𝜕Ω_𝑚/Γ and 𝑢_𝑐=0 on the boundary 𝜕Ω_𝑐/Γ.

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 (⋅,⋅)_𝐷 denote the 𝐿² inner product on the domain 𝐷 (𝐷=Ω_𝑐 or Ω_𝑚) and let ⟨⋅,⋅⟩ denote the 𝐿² inner product on the interface Γ or the duality pairing between (𝐻^1/2₀₀(Γ))^′ and 𝐻^1/2₀₀(Γ) [5]. Define the spaces 𝑋_𝑐=𝑣∈𝐻¹Ω_𝑐^𝑑∣𝑣=0on𝜕Ω_𝑐Γ,𝑄_𝑐=𝐿²Ω_𝑐,𝑋_𝑚=𝜓∈𝐻¹Ω_𝑚∣𝜓=0on𝜕Ω_𝑚Γ, (8) the bilinear forms 𝑎_𝑚𝜙_𝑚,𝜓=𝕂∇𝜙_𝑚,∇𝜓_Ω𝑚,𝑎_𝑐𝑢_𝑐,𝑣=2𝜈𝔻𝑢_𝑐,𝔻(𝑣)_Ω𝑐,𝑏_𝑐𝑣,𝑞=−∇⋅𝑣,𝑞_Ω𝑐, (9) and the projection onto the tangent space on Γ: 𝑃_𝜏𝑢=^𝑑−1_𝑗=1𝑢⋅𝝉_𝑗𝝉_𝑗. (10)

With these notations, the weak formulation of the coupled Navier-Stokes-Darcy model with BJS interface condition is given as follows [46–51]: find (𝑢_𝑐,𝑝_𝑐,𝜙_𝑚)∈𝑋_𝑐×𝑄_𝑐×𝑋_𝑚 such that 𝑢_𝑐⋅∇𝑢_𝑐,𝑣_Ω𝑐+𝑎_𝑐𝑢_𝑐,𝑣+𝑏_𝑐𝑣,𝑝_𝑐−𝑏_𝑐𝑢_𝑐,𝑞+𝑎_𝑚𝜙_𝑚+,𝜓𝑔𝜙_𝑚,𝑣⋅𝑛_𝑐−𝑢_𝑐⋅𝑛_𝑐+√,𝜓𝛼𝜈𝐝(∏)𝑃trace^𝜏𝑢_𝑐,𝑃_𝜏𝑣=𝑓_𝑚,𝜓_Ω𝑚+𝑓_𝑐,𝑣_Ω𝑐+𝑔𝑧,𝑣⋅𝑛_𝑐,∀𝑣,𝑞,𝜓∈𝑋_𝑐×𝑄_𝑐×𝑋_𝑚. (11)

Assume that we have in hand regular subdivisions of Ω_𝑚 and Ω_𝑐 into finite elements with mesh size ℎ. Then one can define finite element spaces 𝑋_𝑚ℎ⊂𝑋_𝑚, 𝑋_𝑐ℎ⊂𝑋_𝑐 and 𝑄_𝑐ℎ⊂𝑄_𝑐. We assume that 𝑋_𝑐ℎ and 𝑄_𝑐ℎ satisfy the inf-sup condition [56, 57] inf_{0≠𝑞ℎ∈𝑄𝑐ℎ}sup_{𝑣0≠ℎ∈𝑋𝑐ℎ}𝑏_𝑐𝑣_ℎ,𝑞_ℎ𝑣_ℎ1𝑞_ℎ0>𝛾, (12) where 𝛾>0 is a constant independent of ℎ. 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 𝑋_𝑚ℎ,  𝑋_𝑐ℎ, and 𝑄_𝑐ℎ that satisfy (12). One example is the Taylor-Hood element pair that we use in the numerical experiments; for that pair, 𝑋_𝑐ℎ consists of continuous piecewise quadratic polynomials and 𝑄_𝑐ℎ 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 (𝑢_𝑐,ℎ,𝑝_𝑐,ℎ,𝜙_𝑚,ℎ)∈𝑋_𝑐ℎ×𝑄_𝑐ℎ×𝑋_𝑚ℎ in the following procedure.(1)The initial value 𝑢⁰_𝑐,ℎ is chosen.(2)For 𝑚=0,1,2,…,𝑀, solve 𝑢^𝑚+1_𝑐,ℎ⋅∇𝑢^𝑚_𝑐,ℎ,𝑣_ℎΩ𝑐+𝑢^𝑚_𝑐,ℎ⋅∇𝑢^𝑚+1_𝑐,ℎ,𝑣_ℎΩ𝑐+𝑎_𝑐𝑢^𝑚+1_𝑐,ℎ,𝑣_ℎ+𝑏_𝑐𝑣_ℎ,𝑝^𝑚+1_𝑐,ℎ−𝑏_𝑐𝑢^𝑚+1_𝑐,ℎ,𝑞_ℎ+𝑎_𝑚𝜙^𝑚+1_𝑚,ℎ,𝜓_ℎ+𝑔𝜙^𝑚+1_𝑚,ℎ,𝑣_ℎ⋅𝑛_𝑐−𝑢^𝑚+1_𝑐,ℎ⋅𝑛_𝑐,𝜓_ℎ+𝛼𝜈√𝐝(∏)𝑃trace^𝜏𝑢^𝑚+1_𝑐,ℎ,𝑃_𝜏𝑣_ℎ=𝑢^𝑚_𝑐,ℎ⋅∇𝑢^𝑚_𝑐,ℎ,𝑣_ℎΩ𝑐+𝑓_𝑚,𝜓_ℎΩ𝑚+𝑓_𝑐,𝑣_ℎΩ𝑐+𝑔𝑧,𝑣_ℎ⋅𝑛_𝑐,∀𝑣_ℎ,𝑞_ℎ,𝜓_ℎ∈𝑋_𝑐ℎ×𝑄_𝑐ℎ×𝑋_𝑚ℎ. (13)(3)Set 𝑢_𝑐,ℎ=𝑢^𝑚+1_𝑐,ℎ,  𝑝_𝑐,ℎ=𝑝^𝑚+1_𝑐,ℎ, and 𝜙_𝑚,ℎ=𝜙^𝑀+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 𝛾_𝑝>0 and a given function 𝜂_𝑝 defined on Γ, 𝛾_𝑝𝕂∇𝜙_𝑚⋅𝑛_𝑚+𝑔𝜙_𝑚=𝜂_𝑝,onΓ. (14) Then, the corresponding weak formulation for the Darcy part is given by the following: for 𝜂_𝑝∈𝐿²(Γ), find 𝜙_𝑚∈𝑋_𝑚 such that 𝑎_𝑚𝜙_𝑚+𝑔𝜙,𝜓^𝑚𝛾_𝑝=𝑓,𝜓_𝑚,𝜓_Ω𝑚+𝜂_𝑝𝛾_𝑝,𝜓,∀𝜓∈𝑋_𝑚. (15)

Second, we can propose the following two Robin-type conditions for the Navier-Stokes equations: for a given constant 𝛾_𝑓>0 and given functions 𝜂_𝑓 and 𝜂_𝑓𝜏 defined on Γ, 𝑛_𝑐⋅𝕋𝑢_𝑐,𝑝_𝑐⋅𝑛_𝑐+𝛾_𝑓𝑢_𝑐⋅𝑛_𝑐=𝜂_𝑓,onΓ,−𝑃_𝜏𝕋𝑢_𝑐,𝑝_𝑐⋅𝑛_𝑐=𝛼𝜈√𝐝(∏)𝑃trace^𝜏𝑢_𝑐,onΓ. (16)

Then, the corresponding weak formulation for the Navier-Stokes equation is given by the following: for 𝜂_𝑓∈𝐿²(Γ), find (𝑢_𝑐,𝑝_𝑐)∈𝑋_𝑐×𝑄_𝑐 such that 𝑢_𝑐⋅∇𝑢_𝑐,𝑣_Ω𝑐+𝑎_𝑐𝑢_𝑐,𝑣+𝑏_𝑐𝑣,𝑝_𝑐−𝑏_𝑐𝑢_𝑐,𝑞+𝛾_𝑓𝑢_𝑐⋅𝑛_𝑐,𝑣⋅𝑢_𝑐+𝛼𝜈√𝐝(∏)𝑃trace^𝜏𝑢_𝑐,𝑃_𝜏𝑣=𝑓_𝑐,𝑣_Ω𝑐+𝜂_𝑓,𝑣⋅𝑛_𝑐,∀𝑣,𝑞∈𝑋_𝑐×𝑄_𝑐. (17)

Our next step is to show that, for appropriate choices of 𝛾_𝑓,  𝛾_p,  𝜂_𝑓, and 𝜂_𝑝, 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 (𝜙_𝑚,𝑢_𝑐,𝑝_𝑐) be the solution of the coupled Navier-Stokes-Darcy system (11) and let (𝜙_𝑚,𝑢_𝑐,𝑝_𝑐) be the solution of the decoupled Navier-Stokes and Darcy equations (15) and (17) with Robin boundary conditions at the interface. Then, (𝜙_𝑚,𝑢_𝑐,𝑝_𝑐)=(𝜙_𝑚,𝑢_𝑐,𝑝_𝑐) if and only if 𝛾_𝑓,  𝛾_𝑝,  𝜂_𝑓,  𝜂_𝑓𝜏, and 𝜂_𝑝 satisfy the following compatibility conditions: 𝜂_𝑝=𝛾_𝑝𝑢_𝑐⋅𝑛_𝑐+𝑔𝜙_𝑚,𝜂_𝑓=𝛾_𝑓𝑢_𝑐⋅𝑛_𝑐−𝑔𝜙_𝑚+𝑔𝑧. (18)(19)

Proof. Adding (15) and (17) together, we obtain the following: given 𝜂_𝑝,𝜂_𝑓∈𝐿²(Γ), find (𝜙_𝑚,𝑢_𝑓,𝑝_𝑐)∈𝑋_m×𝑋_𝑐×𝑄_𝑐 such that 𝑢_𝑐⋅∇𝑢_𝑐,𝑣_Ω𝑐+𝑎_𝑐𝑢_𝑐,𝑣+𝑏_𝑐𝑣,𝑝_𝑐−𝑏_𝑐𝑢_𝑐,𝑞+𝑎_𝑚𝜙_𝑚,𝜓+𝛾_𝑓𝑢_𝑐⋅𝑛_𝑐,𝑣⋅𝑛_𝑐+𝑔𝜙_𝑚𝛾_𝑝+√,𝜓𝛼𝜈𝐝(∏)𝑃trace^𝜏𝑢_𝑐,𝑃_𝜏𝑣=𝑓_𝑚,𝜓_Ω𝑚+𝑓_𝑐,𝑣_Ω𝑐+𝜂_𝑓,𝑣⋅𝑛_𝑐+𝜂_𝑝𝛾_𝑝,𝜓,∀𝑣,𝑞,𝜓∈𝑋_𝑚×𝑋_𝑐×𝑄_𝑐. (20)
For the necessity of the lemma, we pick 𝜓=0 and 𝑣 such that 𝑃_𝜏𝑣=0 in (11) and (20); then by subtracting (20) from (11), we get 𝜂_𝑓−𝛾_𝑓𝑣_𝑓⋅𝑛_𝑐+𝑔𝜙_𝑚−𝑔𝑧,𝑣⋅𝑛_𝑐∀=0,𝑣∈𝑋_𝑐with𝑃_𝜏𝑣=0, (21) 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 (𝜙_𝑚,𝑢_𝑐,𝑝_𝑐) 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 𝜂⁰_𝑝 and 𝜂⁰_𝑓 are guessed. They may be taken to be zero.(2) For 𝑘=0,1,2,…, independently solve the Darcy and Navier-Stokes equations constructed above. More precisely, 𝜙^𝑘_𝑚∈𝑋_𝑚 is computed from 𝑎_𝑚𝜙^𝑘_𝑚,𝜓+𝑔𝜙^𝑘_𝑚𝛾_𝑝=𝜂,𝜓^𝑘𝑝𝛾_𝑝,𝜓+𝑓_𝑚,𝜓_Ω𝑚,∀𝜓∈𝑋_𝑚, (22) and 𝑢^𝑘_𝑐∈𝑋_𝑐 and 𝑝^𝑘_𝑐∈𝑄_𝑐 are computed from the following Newton iteration.(i) Initial value 𝑢_𝑐^𝑘,0 is chosen for the Newton iteration. For instance, it may be taken to be 𝑢_𝑐^0,0=0 and 𝑢_𝑐^𝑘,0=𝑢_𝑐^𝑘−1 for 𝑘=1,2,….(ii) For 𝑚=0,1,2,…,𝑀, solve 𝑢_𝑐^𝑘,𝑚+1⋅∇𝑢_𝑐^𝑘,𝑚,𝑣_Ω𝑐+𝑢_𝑐^𝑘,𝑚⋅∇𝑢_𝑐^𝑘,𝑚+1,𝑣_Ω𝑐+𝑎_𝑐𝑢_𝑐^𝑘,𝑚+1,𝑣+𝑏_𝑐𝑣,𝑝_𝑐^𝑘,𝑚+1−𝑏_𝑐𝑢_𝑐^𝑘,𝑚+1,𝑞+𝛾_𝑓𝑢_𝑐^𝑘,𝑚+1⋅𝑛_𝑐,𝑣⋅𝑛_𝑐+𝛼𝜈√𝐝(∏)𝑃trace^𝜏𝑢_𝑐^𝑘,𝑚+1,𝑃_𝜏𝑣=𝑢_𝑐^𝑘,𝑚⋅∇𝑢_𝑐^𝑘,𝑚,𝑣_Ω𝑐+𝜂^𝑘_𝑓,𝑣⋅𝑛_𝑐+𝑓_𝑐,𝑣_Ω𝑐,∀𝑣,𝑞,𝜓∈𝑋_𝑐×𝑄_𝑐×𝑋_𝑚. (23)(iii) Set 𝑢^𝑘_𝑐=𝑢_𝑐^𝑘,𝑀+1 and 𝑝^𝑘_𝑐=𝑝_𝑐^𝑘,𝑀+1.(3)𝜂_𝑝^𝑘+1 and 𝜂_𝑓^𝑘+1 are updated in the following manner: 𝜂_𝑓^𝑘+1=𝛾_𝑓𝛾_𝑝𝜂^𝑘_𝑝−𝛾1+^𝑓𝛾_𝑝𝑔𝜙^𝑘_𝑚𝜂+𝑔𝑧,_𝑝𝑘+1=−𝜂^𝑘_𝑓+𝛾_𝑓+𝛾_𝑝𝑢^𝑘_𝑐⋅𝑛_𝑐+𝑔𝑧. (24)

Then the corresponding domain decomposition finite element method is proposed as follows.(1) Initial values 𝜂⁰_𝑝,ℎ and 𝜂⁰_𝑓,ℎ are guessed. They may be taken to be zero. (2) For 𝑘=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, 𝜙^𝑘_𝑚,ℎ∈𝑋_𝑚ℎ is computed from 𝑎_𝑚𝜙^𝑘_𝑚,ℎ,𝜓_ℎ+𝑔𝜙^𝑘_𝑚,ℎ𝛾_𝑝,𝜓_ℎ=𝜂^𝑘_𝑝,ℎ𝛾_𝑝,𝜓_ℎ+𝑓_𝑚,𝜓_ℎΩ𝑚,∀𝜓_ℎ∈𝑋_𝑚ℎ, (25) and 𝑢^𝑘_𝑐,ℎ∈𝑋_𝑐ℎ and 𝑝^𝑘_𝑐,ℎ∈𝑄_𝑐ℎ are computed from the following Newton iteration.(i) Initial value 𝑢^𝑘,0_𝑐,ℎ is chosen for the Newton iteration. For instance, it may be taken to be 𝑢^0,0_𝑐,ℎ=0 and 𝑢^𝑘,0_𝑐,ℎ=𝑢^𝑘−1_𝑐,ℎ for 𝑘=1,2,….(ii) For 𝑚=0,1,2,…,𝑀, solve 𝑢^𝑘,𝑚+1_𝑐,ℎ⋅∇𝑢^𝑘,𝑚_𝑐,ℎ,𝑣_ℎΩ𝑐+𝑢^𝑘,𝑚_𝑐,ℎ⋅∇𝑢^𝑘,𝑚+1_𝑐,ℎ,𝑣_ℎΩ𝑐+𝑎_𝑐𝑢^𝑘,𝑚+1_𝑐,ℎ,𝑣_ℎ+𝑏_𝑐𝑣_ℎ,𝑝_𝑐^𝑘,𝑚+1−𝑏_𝑐𝑢^𝑘,𝑚+1_𝑐,ℎ,𝑞_ℎ+𝛾_𝑓𝑢^𝑘,𝑚+1_𝑐,ℎ⋅𝑛_𝑐,𝑣_ℎ⋅𝑛_𝑐+𝛼𝜈√𝐝(∏)𝑃trace^𝜏𝑢^𝑘,𝑚+1_𝑐,ℎ,𝑃_𝜏𝑣_ℎ=𝑢^𝑘,𝑚_𝑐,ℎ⋅∇𝑢^𝑘,𝑚_𝑐,ℎ,𝑣_ℎΩ𝑐+𝜂^𝑘_𝑓,ℎ,𝑣_ℎ⋅𝑛_𝑐+𝑓_𝑐,𝑣_ℎΩ𝑐,∀𝑣_ℎ,𝑞_ℎ,𝜓_ℎ∈𝑋_𝑐ℎ×𝑄_𝑐ℎ×𝑋_𝑚ℎ. (26)(iii) Set 𝑢^𝑘_𝑐,ℎ=𝑢^𝑘,𝑚+1_𝑐,ℎ  and  𝑝^𝑘_𝑐,ℎ=𝑝^𝑘,𝑀+1_𝑐,ℎ.(3)𝜂^𝑘+1_𝑝,ℎ and 𝜂^𝑘+1_𝑓,ℎ are updated in the following manner: 𝜂^𝑘+1_𝑓,ℎ=𝛾_𝑓𝛾_𝑝𝜂^𝑘_𝑝,ℎ−𝛾1+^𝑓𝛾_𝑝𝑔𝜙^𝑘_𝑚,ℎ𝜂+𝑔𝑧,_{𝑘+1𝑝,ℎ}=−𝜂^𝑘_𝑓,ℎ+𝛾_𝑓+𝛾_𝑝𝑢^𝑘_𝑐,ℎ⋅𝑛_𝑐+𝑔𝑧. (27)

5. Numerical Example

Example 1. Consider the model problem (2)–(6) with the BJS interface condition (7) on Ω=[0,𝜋]×[−1,1] with Ω_𝑚=[0,𝜋]×[0,1] and Ω_𝑐=[0,𝜋]×[−1,0]. Choose (𝛼𝜈√𝐝/trace(∏))=1, 𝜈=1, 𝑔=1, 𝑧=0, and 𝕂=𝐾𝕀, where 𝕀 is the identity matrix and 𝐾=1. The boundary condition data functions and the source terms are chosen such that the exact solution is given by 𝜙_𝑚=𝑒^𝑦−𝑒^−𝑦sin(𝑥)𝑒^𝑡,𝑢_𝑐=𝐾𝜋sin2𝜋𝑦cos(𝑥)𝑒^𝑡,𝐾−2𝐾+𝜋²sin²𝜋𝑦sin(𝑥)𝑒^𝑡𝑇,𝑝_𝑐=0. (28) We divide Ω_𝑚 and Ω_𝑐 into rectangles of height ℎ=1/𝑁 and width 𝜋ℎ, where 𝑁 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, Table 1 provides errors for different choices of ℎ. Using linear regression, the errors in Table 1 satisfy 𝑢_𝑐,ℎ−𝑢_𝑐0≈0.714ℎ^3.011,𝑢_𝑐,ℎ−𝑢_𝑐1≈3.867ℎ^1.987,𝑝_𝑐,ℎ−𝑝_𝑐0≈5.123ℎ^3.129,𝜙_𝑚,ℎ−𝜙_𝑚0≈0.354ℎ^2.998,𝜙_𝑚,ℎ−𝜙_𝑚1≈1.556ℎ^1.995. (29) 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 𝐿² norm of 𝑢_𝑐 and 𝜙_𝑚, second order with respect to 𝐻¹ norm of 𝑢_𝑐 and 𝜙_𝑚, and second order with respect to 𝐿² norms of 𝑝_𝑐. In particular, the third-order convergence rate of 𝑝_𝑐 observed above, which is better than the approximation capability of the linear element, is mainly due to the special analytic solution 𝑝=0.

Table 1 

Errors of the finite element method for the steady Navier-Stokes-Darcy model with BJS interface condition.

For the parallel DDM with 𝜈=1,  𝐾=1,  𝛾_𝑓=0.3, and ℎ=1/32, Figures 2 and 3 show the 𝐿² errors of hydraulic head, velocity, pressure, and 𝜂_𝑓. We can see that the parallel domain decomposition method is convergent for 𝛾_𝑓≤𝛾_𝑝. Moreover, Figures 4 and 5 show that a smaller 𝛾_𝑓/𝛾_𝑝 leads to faster convergence.

(a)

(b)

Figure 2 

Convergence for the velocity of the free flow (a) and the hydraulic head of the porous medium flow (b) versus the iteration counter 𝑚 for the parallel DDM with BJS interface condition.

(a)

(b)

Figure 3 

Convergence for the pressure of the free flow (a) and 𝜂𝑓 (b) versus the iteration counter 𝑚 for the parallel DDM with BJS interface condition.

(a)

(b)

Figure 4 

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)

Figure 5 

Geometric convergence rate of the pressure of the free flow (a) and 𝜂𝑓 (b) versus the iteration counter 𝑚 for the parallel DDM with BJS interface condition.

Then Tables 2 and 3 list some 𝐿² errors in velocity, hydraulic head, pressure, and 𝜂_𝑓 for the parallel domain decomposition method with 𝛾_𝑓=0.3 and 𝛾_𝑝=1.2. The data in these two tables indicate the geometric convergence rate 𝛾_𝑓/𝛾_𝑝 since all the error ratios are less than (𝛾_𝑓/𝛾_𝑝)⁴=(1/4)⁴=0.0625.

Table 2 

𝐿² errors in velocity and hydraulic head for the parallel DDM with BJS interface condition.

Table 3 

𝐿² errors in pressure and 𝜂_𝑓 for the parallel DDM with BJS interface condition.

Finally, for the preconditioning feature of the domain decomposition method, Table 4 shows the number of iterations 𝑀 is independent of the grid size ℎ. Here, we set 𝛾_𝑆=0.3,  𝛾_𝐷=1.2,  𝜈=1, and 𝐾=1. Let 𝜙^𝑘_ℎ,  𝑢^𝑘_ℎ, and 𝑝^𝑘_ℎ denote the finite element solutions of 𝜙^𝑘_𝐷,  𝑢^𝑘_𝑆, and 𝑝^𝑘_𝑆 at the 𝑘th step of the domain decomposition algorithm. The criterion used to stop the iteration, that is, to determine the value 𝑀, is ‖𝑢^𝑘_ℎ−𝑢_ℎ^𝑘−1‖₀+‖𝜙^𝑘_ℎ−𝜙_ℎ^𝑘−1‖₀+‖𝑝^𝑘_ℎ−𝑝_ℎ^𝑘−1‖₀<𝜀, where the tolerance 𝜀=10⁻⁵.

Table 4 

The iteration counter 𝑀 versus the grid size ℎ for both the parallel Robin-Robin domain decomposition method with BJS interface condition.

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.

Acknowledgments

This work is partially supported by DOE Grant DE-FE0009843, National Natural Science Foundation of China (11175052).