# Analysis of the diffuse interface method for the Stokes-Darcy coupled problem

We consider the interaction between a free flowing fluid and a porous medium flow, where the free flowing fluid is described using the time dependent Stokes equations, and the porous medium flow is described using Darcy's law in the primal formulation. To solve this problem numerically, we use the diffuse interface approach, where the weak form of the coupled problem is written on an extended domain which contains both Stokes and Darcy regions. This is achieved using a phase-field function which equals one in the Stokes region and zero in the Darcy region, and smoothly transitions between these two values on a diffuse region of width ϵ around the Stokes-Darcy interface. We prove the convergence of the diffuse interface formulation to the standard, sharp interface formulation, and derive the rates of the convergence. This is performed by analyzing the modeling error of the diffuse interface approach at the continuous level, and by deriving the a priori error estimates for the diffuse interface method at the discrete level. The convergence rates are also derived computationally in a numerical example.

## Authors

• 3 publications
• 1 publication
04/06/2021

### Analysis of the Stokes-Darcy problem with generalised interface conditions

Fluid flows in coupled systems consisting of a free-flow region and the ...
06/19/2021

### Sharp-interface limits of the diffuse interface model for two-phase inductionless magnetohydrodynamic fluids

In this paper, we propose and analyze a diffuse interface model for indu...
07/03/2020

### A global-in-time domain decomposition method for the coupled nonlinear Stokes and Darcy flows

We study a decoupling iterative algorithm based on domain decomposition ...
11/18/2021

### Parabolic interface reconstruction for 2D volume of fluid methods

For capillary driven flow the interface curvature is essential in the mo...
12/24/2021

### Computing Viscous Flow Along a 3D Open Tube Using the Immerse Interface Method

In a companion study <cit.>, we present a numerical method for simulatin...
01/15/2020

### Robust preconditioning for coupled Stokes-Darcy problems with the Darcy problem in primal form

The coupled Darcy-Stokes problem is widely used for modeling fluid trans...
09/01/2020

### A hybrid WENO method with modified ghost fluid method for compressible two-medium flow problems

In this paper, we develop a simplified hybrid weighted essentially non-o...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The interaction between a free flowing fluid and a porous medium flow, commonly formulated as a Stokes-Darcy coupled system, has been used to describe problems arising in many applications, including environmental sciences, hydrology, petroleum engineering and biomedical engineering. Hence, the development of numerical methods for Stokes-Darcy problems has been an active area of research. Most existing numerical methods are based on a sharp interface approach, in the sense that the interface between the Stokes and Darcy regions is parametrized using an exact specification of its geometry and location, and the nodes in the computational mesh align with the interface. They include both monolithic and partitioned numerical methods. Some of the recently developed monolithic schemes include a two-grid method with backtracking for the stationary monolithic Stokes-Darcy problem proposed in [1] and a mortar multiscale finite element method presented in [2]. We also mention the work in [3, 4] based on the Nitsche’s penalty method. Non-iterative partitioned schemes based on various time-discretization strategies were developed in [5, 6, 7, 8, 9, 10]. A third-order in time implicit-explicit algorithm based on the combination of the Adams–Moulton and the Adams–Bashforth scheme was proposed in [11], and iterative domain decomposition methods based on generalized Robin coupling conditions were derived in [12, 13].

While the sharp interface methods are widely used, the explicit interface parametrization may be difficult to obtain in case of complex geometries. The exact location is sometimes not known exactly or the geometry is complicated, making a proper approximation of the integrals error-prone and difficult to automate. This often occurs when geometries are obtained implicitly using imaging data, commonly used in patient-specific biomedical simulations. Hence, as one alternative to sharp interface approaches, diffuse interface methods have been introduced [14, 15, 16, 17, 18, 19, 20]. They are also known as the phase-field or the diffuse domain methods. While they have features in common with the level set method [21, 22, 23, 24], they are fundamentally different since the level set method tracks the exact sharp interface without introducing any diffuse layers [20]. Other conceptually similar approaches include the fictitious domain method with a spread interface [25, 26], the immersed boundary method with interface forcing functions based on Dirac distributions [27, 28] and the fat boundary method [29, 30], among others. The diffuse interface method received strong attention from the applications point of view [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41], however, many difficult questions remain to be theoretically addressed both at the continuous and the discrete level. One of the fundamental mathematical questions, whether the diffuse interface converges to a sharp interface when the width of the diffuse interface tends to zero, remains an open problem for many phase-field models.

Theoretically, the diffuse interface approach has been widely studied for elliptic problems and two-phase flow problems. Analysis of the diffuse interface method for elliptic problems has been studied in [18, 42, 43, 44, 16]. An elliptic problem was also considered in [45], where diffuse formulations of Nitsche’s method for imposing Dirichlet boundary conditions on phase-field approximations of sharp domains were studied. A Cahn-Larche phase-field model approximating an elasticity sharp interface problem has been considered in [46]. Advection–diffusion on evolving diffuse interfaces has been studied in [33, 32]. We also mention the work in [47] where the interaction between the two-phase free flow and two-phase porous media flow was considered, but the diffuse interface method was used to separate different phases in each region, while the Stokes-Darcy interface was captured by a mesh aligned with the interface and the two problems were decoupled using a classical approach. Diffuse interface models for two-phase flows have been extensively analyzed, see, e.g., [48, 49, 50, 51, 52, 53, 54]. However, only recently a rigorous sharp interface limit convergence result with convergence rates in strong norms has been proved [55]. The main challenges in the analysis of the diffuse interface problems involve proving the convergence of the diffuse interface to the sharp interface, and the estimation of the modeling error, which can be quite difficult to analyze and remains open for many problems

We consider the fluid–porous medium interaction described by the time–dependent Stokes–Darcy coupled problem. We discretize the problem in time using the backward Euler method, and in space using the finite element method. The focus of this paper is the convergence analysis of the diffuse interface to the sharp interface as the width of the diffuse layer goes to zero. This is performed in two steps. First, we derive the error estimates and prove the convergence for the finite element approximation of the diffuse interface method. Then, we analyze the convergence of the continuous diffuse interface formulation to the continuous sharp interface formulation, obtaining the convergence rates with respect to the width of the diffuse layer. The rates of convergence are also computed numerically.

The rest of this paper is organized as follows. Section 2 presents the mathematical model and introduces a diffuse interface formulation. We also prove that the problem is well-posed, and discretize the problem in time and space. The error analysis of the discrete diffuse interface problem is performed in Section 3, and the analysis of the modeling error is performed in Section 4. Numerical examples are presented in Section 5, and conclusions are drawn in Section 6.

## 2 The mathematical model

### 2.1 The sharp interface model

Let denote the fluid domain and denote the porous medium domain. We assume that are open, smooth sets of the same dimension, and that the fluid–porous medium interface is the common boundary between the two domains, i.e. To model the fluid flow, we use the time–dependent Stokes equations, given as follows:

 ρ∂tu=∇⋅σ(∇u,π)+FinΩF, (1) ∇⋅u=0inΩF, (2)

where is the fluid velocity, is the fluid density,

is the fluid Cauchy stress tensor,

is the fluid viscosity, denotes the pressure and is the density of volume force.

The porous medium flow is governed by Darcy’s law, written in the primal formulation as:

 c0∂tp−∇⋅(κ∇p)=ginΩD, (3)

where is the Darcy pressure, is the storativity and is the hydraulic conductivity tensor. We assume that is uniformly bounded and positive definite, so that where , and is the spectrum of . We note that the porous medium flow can be computed from as

Coupling conditions: The Stokes and Darcy problems are coupled using the following interface conditions:

• The conservation of the mass:

 u⋅n=−κ∇p⋅nonΓ×(0,T),

where is the unit normal to the fluid domain boundary.

• The Beavers-Joseph-Saffman condition:

 αBJu⋅τi+σ(u,π)n⋅τi=0,i=1,2,…,d−1,onΓ×(0,T),

where , , is the orthonormal basis for the tangential space of and is the Beavers-Joseph-Saffman-Jones coefficient.

• The balance of pressure:

 −σ(u,π)n⋅n=ponΓ×(0,T).

Boundary conditions: We split the boundaries as: , , and prescribe the following boundary conditions:

 u=0onΓ1F, σn=0onΓ2F, p=0onΓ1D, ∂np=0onΓ2D.
###### Remark

In this paper, for simplicity, we consider only homogeneous boundary conditions. However, inhomogeneous boundary conditions can be easily included in the analysis by using the appropriate extension operators to homogenize the boundary conditions.

#### 2.1.1 Weak formulation

Given an open set , we consider the usual Sobolev spaces , with We introduce the following functional spaces:

 V={u∈(H1ΩF))d:u=0onΓ1F}, X={ψ∈H1(ΩD):ψ=0onΓ1D}.

The weak solution to the sharp interface Stokes-Darcy problem is defined as follows.

###### Definition 1 (Weak solution of the sharp interface form)

We say that is a weak solution if, for every , the following equality is satisfied in :

 ddt(∫ΩFρu⋅v+c0∫ΩDpψ)+2ν∫ΩFD(u):D(v)+∫ΩDκ∇p⋅∇ψ−∫ΩF(∇⋅v)π+∫ΩF(∇⋅u)φ −∫Γψu⋅n+∫Γpv⋅n+αBJd−1∑i=1∫Γ(u⋅τi)(v⋅τi)=∫ΩFF⋅v+∫ΩDgψ.

Here, denotes the standard space of distributions on . The weak solution satisfies the following energy equality:

 12ddt(ρ∥u∥2L2(ΩF)+c0∥p∥2L2(ΩD))+2ν∥D(u)∥2L2(ΩF)+αBJd−1∑i=1∥u⋅τi∥2L2(Γ)+∥κ12∇p∥2L2(ΩD)=0.

The theory of weak solutions for the Stokes-Darcy system is by now well understood, see, e.g., [56] for the stationary case and [57] for the evolution case.

### 2.2 The diffuse interface formulation

Denote by the Heaviside function which equals one in and zero in . Let a phase field function be a regularization of the Heaviside function such that in and in , and transitions between these two values on a “diffuse” layer of width . We suppose that and . More precisely, we assume that the phase field function satisfies conditions from [18] (see Section 3, conditions (S1)-(S3)). An example of such function is , where is a signed distance function, and

 S(t)=⎧⎪⎨⎪⎩t,t≤1,t|t|,t>1.

Then, we can write:

 ∫ΩFFdΩ=∫ΩFχdΩ≈∫ΩFΦϵdΩ, ∫ΩDFdΩ=∫ΩF(1−χ)dΩ≈∫ΩF(1−Φϵ)dΩ, ∫ΓfdΓ=∫ΩfδΓdΩ≈∫Ωf|∇Φϵ|dΩ,

where is a Dirac distribution at the interface . Using the approximations above, and , we can approximate the interface integrals as:

 ∫Γψu⋅n≈−∫Ωψu⋅∇Φϵ,∫Γpv⋅n≈−∫Ωpv⋅∇Φϵ

and

where the tangent vectors

are obtained directly from the phase-field function using the algorithm described in [58, 31]. We introduce the following regularization of

 Φϵ=(1−2δ)Φϵ+δ, (4)

where is a small positive number. Therefore, and . Since we will keep fixed throughout the next two sections, we do not write explicitly dependence of on .

Every weight, , induces a measure with density , over the Borel sets of . For simplicity, this measure will also be denoted by . For a Borel set , we let . Similar as in [18], we define the weighted space associated to the measure by:

 Lp(Ω,ω)={ψ:Ω→R:|ψ|pω∈L1(Ω)},

with norm

 ∥ψ∥pLp(Ω,ω)=∫Ω|ψ|pωdx.

Associated with weighted spaces, we define the weighted Sobolev spaces:

 Hk(Ω,ω)={ψ∈L2(Ω,ω):Dαψ∈L2(Ω,ω),∀|α|≤k),

with the corresponding norm

 ∥ϕ∥2Hk(Ω,ω)=∑|α|≤k∥Dαψ∥2L2(Ω,ω).

Notice that due to the regularization (4), the weighted spaces with weight are the same as the classical spaces as shown in the following lemma:

###### Lemma 1

For , we have , with equivalent norms given by:

 √δ∥f∥L2(Ω)≤∥f∥L2(Ω,Φϵ)≤∥f∥L2(Ω).

Similarly, , with equivalent norms given by:

 √δ∥f∥H1(Ω)≤∥f∥H1(Ω,Φϵ)≤∥f∥H1(Ω).

###### Proof

Notice that, due to the regularization (4), the phase field function is bounded as . Therefore, we have:

 ∥f∥2L2(Ω)=∫Ωf2≤1δ∫ΩΦϵf2=1δ∥f∥2L2(Ω,Φϵ), ∥f∥2L2(Ω,Φϵ)=∫ΩΦϵf2≤∫Ωf2=∥f∥2L2(Ω).

Because in the phase field formulation all unknowns are defined on the whole domain, , we define the following Hilbert spaces for :

 Vk,ϵ={u∈(Hk(Ω,Φϵ))d:u=0onΓ1F},Vϵ=V1,ϵ. Xk,ϵ={ψ∈Hk(Ω,1−Φϵ):ψ=0onΓ1D},Xϵ=X1,ϵ.

Space is associated with the fluid velocity, while space is associated with the porous pressure. The following results from [18] and [59] will be used.

###### Lemma 2 (Trace inequality)

Let be sufficiently small. Then, there exists a constant such that for and for , we have:

 ∫Ω|v|p|∇Φϵ|≤CT∥v∥pW1,p(Ω,Φϵ).

###### Lemma 3 (Poincaré inequality)

We define as the closure of in and set On these spaces, the following Poincaré inequality holds:

 ∥v∥Lp(Ω,Φϵ)≤CP∥∇v∥Lp(Ω,Φϵ),∀v∈W1,p0(Ω,Φϵ),

where depends on the diameter of , but is independent of .

###### Definition 2 (Weak solution of the diffuse interface form)

We say that is a weak solution if, for every , the following equality is satisfied in :

 −∫Ω(∇⋅vϵ)πΦϵ+∫Ω(∇⋅uϵ)φΦϵ+∫Ωψuϵ⋅∇Φϵ−∫Ωpϵv⋅∇Φϵ (5) +αBJd−1∑i=1∫Ω(u⋅~τi)(v⋅~τi)|∇Φϵ|=∫ΩF⋅vΦϵ+∫Ωgψ(1−Φϵ).

The goal of this paper is to analyze the convergence of the discrete phase field formulation of the Stokes-Darcy system to the solution of the continuous sharp interface problem. This will be done in two steps. In the first step, in Section 3, we fix and the phase field function, and prove the error estimates for the finite element approximation of the diffuse interface formulation (5). In the second step, in Section 4, we analyze the convergence of the continuous diffuse interface formulation to the continuous sharp interface formulation.

Since in this section and Section 3 parameter is fixed, we will omit writing superscript in the unknowns for simplicity of notation. We start by proving the well-posedness for the phase field formulation.

###### Theorem 1

There exists a unique weak solution in sense of Definition 2 to the Stokes-Darcy phase field formulation problem.

###### Proof

Because of the regularization (4) the proof is similar to the well-posedness proof for the Stokes-Darcy system (e.g. [57, 60]). Therefore, here we just outline the main steps of the proof. We define a function space

 Vϵdiv:={u∈Vϵ:∇⋅u=0}.

Notice that the solution is an element of . Now, we define a bilinear form on as:

 a((u,p),(v,ψ)):=2ν∫ΩD(u):D(v)Φ+∫Ωκ∇p⋅∇ψ(1−Φ)

and a linear functional as:

 ⟨F,(v,ψ)⟩:=∫ΩF⋅vΦ+∫Ωgψ(1−Φ).

First we prove the following lemma.

###### Lemma 4

The bilinear form is continuous and coercive on .

###### Proof

In order to show the continuity, special attention has to be given to the terms arising from the coupling at the interface. In particular, using the Cauchy-Schwarz inequality followed by Lemma 2, we obtain:

 ∫Ωpv⋅∇Φ ≤(∫Ω|p|2|∇Φ|)12(∫Ω|v|2|∇Φ|)12≤(∫Ω|p|2|∇(1−Φ)|)12(∫Ω|v|2|∇Φ|)12 ≤CT∥p∥H1(Ω,1−Φ)∥v∥H1(Ω,Φ).

Similarly, we have

 ∫Ωψu⋅∇Φ≤CT∥ψ∥H1(Ω,1−Φ)∥u∥H1(Ω,Φ),

and

 d−1∑i=1∫Ω(u⋅~τi)(v⋅~τi)|∇Φ| ≤d−1∑i=1CT∥u⋅~τi∥H1(Ω,Φ)∥v⋅~τi∥H1(Ω,Φ)≤CT∥u∥H1(Ω,Φ)∥v∥H1(Ω,Φ).

Other terms can be bounded in the standard way. Therefore, we obtain the following estimate:

 a((u,p),(v,ψ))≤ 2ν∥D(u)∥L2(Ω,Φ)∥D(v)∥L2(Ω,Φ)+∥κ∇p∥L2(Ω,1−Φ)∥∇ψ∥L2(Ω,1−Φ) +CT∥p∥H1(Ω,1−Φ)∥v∥H1(Ω,Φ)+CT∥ψ∥H1(Ω,1−Φ)∥u∥H1(Ω,Φ) +αBJCT∥u∥2H1(Ω,Φ)∥v∥2H1(Ω,Φ) ≤ D∥(u,p)∥Vϵdiv×Xϵ∥(v,ψ)∥Vϵdiv×Xϵ

where

 D=max{αBJCT,2ν,k∗,CT}.

To show coercivity, we take to obtain:

 a((u,p),(u,p))= 2ν∫ΩD(u):D(u)Φ+∫Ωκ∇p⋅∇p(1−Φ)+αBJd−1∑i=1∫Ω(u⋅~τi)2|∇Φ| ≥ 2ν∥D(u)∥2L2(Ω,Φ)+∥κ12∇p∥2L2(Ω,1−Φ).

Using the standard Korn and Poincare inequalities on together with Lemma 1, we have

 a((u,p),(u,p))≥D∥(u,p),(u,p)∥2Vϵdiv×Xϵ,

where

 D=min{νδC2K,νδC2PC2K,k∗2,k∗δ2C2PC2K},

which completes the proof of Lemma 4.

From estimates analogous to the ones used in Lemma 4, it follows that . Therefore, by using standard methods (e.g. Galerkin method), one can easily prove that there exists a unique solution to the following evolution problem in :

 ddt(∫Ωρu⋅vΦ+c0∫Ωpψ(1−Φ))+a((u,p),(v,ψ))=⟨F,(v,ψ)⟩,(v,ψ)∈Vϵdiv×Xϵ. (6)

To finish the proof it remains to construct the pressure corresponding to the solution of (6). This follows in analogous way as in the standard theory of the Navier-Stokes equation. We observe that is a surjective operator. Namely, because of the regularization (4), the equation is equivalent to . Therefore, the surjectivity follows directly from the standard Bogovskiǐ theorem (e.g. [61, Section 3.3]). Now, the existence of pressure can be proved in the same way as in the case of the Navier-Stokes equations (see, e.g., [62, Proposition III.1.1.]).

To discretize the problem (5) in space, we use the finite element method. Let be a conforming triangulation of and let , , . We assume that the mesh is regular, that the is a measure of the grid size, and that the fluid and pressure spaces and satisfy the discrete inf-sup condition necessary to ensure the stability of the finite element discretization [63]. We also assume that and include piecewise polynomials of degree at least and that includes piecewise polynomials of degree .

Let for , where denotes the time step, and is the final time. For time discretization, we use the backward Euler method. Let denote the approximation of a time-dependent function at time level . We will use the following notation for the discrete time derivative:

 dtψn+1=ψn+1−ψnΔt.

The discrete problem is given by: Find and such that for every and , we have:

 ρ∫Ωdtun+1h⋅vhΦ+c0∫Ωdtpn+1hψh(1−Φ)+2ν∫ΩD(un+1h):D(vh)Φ+∫Ωκ∇pn+1h⋅∇ψh(1−Φ) −∫Ωπn+1h∇⋅vhΦ+∫Ωφh(∇⋅un+1h)Φ+∫Ωψhun+1h⋅∇Φ−∫Ωpn+1hvh⋅∇Φ (7) +αBJd−1∑i=1∫Ω(un+1h⋅~τi)(vh⋅~τi)|∇Φ|=∫ΩFh⋅vhΦ+∫Ωghψh(1−Φ).

## 3 Error analysis

We assume that for the chosen finite element spaces and , there exists a projection operator satisfying:

 ∫Ω∇⋅(u−Phu)φΦ=0 ∀φ∈XϵF,h,u∈Vϵ, (8) ∥u−Phu∥H1(Ω,Φ)≤Cδ−12hk∥u∥Hk+1(Ω,Φ) ∀u∈Vk+1, (9)

where and C is a positive constant independent of and . Details about the existence of the projection operator in case when can be found in [63] for and in [64] for . Using those results, the relations (8)-(9) follow from Lemma 1. In this section, we assume that is a fixed, positive number.

Let be a projection operator onto such that

 ∥p−Shp∥H1(Ω,1−Φ)≤Cδ−12hk∥p∥Hk+1(Ω,1−Φ)∀p∈Xk+1,ϵD, (10)

and let be a projection operator onto such that

 ∥π−Rhπ∥L2(Ω,Φ)≤Cδ−12hs+1∥π∥Hs+1(Ω,Φ)∀π∈Xk+1,ϵF, (11)

where .

Let denote that there exists a positive constant , independent , and , such that . We introduce the following time discrete norms:

where . Note that they are equivalent to the continuous norms since we use piecewise constant approximations in time. Furthermore, the following inequality holds:

 ΔtN−1∑n=1∥dtvn+1∥2X≲∥∂tv∥2L2(0,T;X).

The main result of this section is stated in the following theorem.

###### Theorem 2

Let be a solution of (5), and let be a solution of (7), with discrete initial data

 (u0h,p0h,π0h)=(Phu0,Shp0,Rhπ0).

Then, the following estimate holds:

 ρ2∥uN−uNh∥2L2(Ω,Φ)+c02∥pN−pNh∥2L2(Ω,1−Φ)+νΔtN−1∑n=0∥D(un+1−un+1h)∥2L2(Ω,Φ) +k∗Δt2N−1∑n=0∥∇(pn+1−pn+1h)∥2L2(Ω,1−Φ) ≲Δt2(ρ2FC2PC2Kνδ∥∂ttu∥2L2(0,T;L2(Ω,Φ))+c20C2Pk∗∥∂ttp∥2L2(0,T;L2(Ω,1−Φ))) +h2k(ρδ∥u∥2L∞(0,T;Hk+1(Ω,Φ))+ρ2C2PC2Kνδ2∥∂tu∥2L2(0,T;Hk+1(Ω,Φ)) +(νδ+C4T(C2P+1)k∗δ+α2BJC4TC2Kδ2(C2P+1))∥u∥2L2(0,T;Hk+1(Ω,Φ))+c0δ∥p∥2L∞(0,T;Hk+1(Ω,1−Φ)) +(c20C2Pk∗δ+(k∗)2k∗δ+C4TC2Kνδ2(C2P+1))∥p∥2L2(0,T;Hk+1(Ω,1−Φ))) +h2s+2dC2Kνδ2∥π∥2L2(0,T;Hs+1(Ω,Φ)). (12)

###### Proof

Define and . We begin by subtracting (7) from (5), obtaining the following error equation:

 ρ∫Ωdten+1U⋅vhΦ+2ν∫ΩD(en+1U):D(vh)Φ−∫Ωen+1π∇⋅vhΦ+∫Ωφh∇⋅en+1UΦ +c0∫Ωdten+1Pψh(1−Φ)+∫Ωκ∇en+1P⋅∇ψh(1−Φ)−∫Ωen+1Pvh⋅∇Φ +∫Ωψhen+1U⋅∇Φ+αBJd−1∑i=1∫Ω(en+1U⋅~τi)(vh⋅~τi)|∇Φ|=R(vhψh), (13)

where

 R(vh,ψh)=ρF∫Ω(dtun+1−∂tun+1)⋅vhΦ+c0∫Ω(dtpn+1−∂tpn+1)ψh(1−Φ).

We decompose the error into the sum of the interpolation error and the approximation error as:

 eπ=π−Rhπ+Rhπ−πh=θπ+δπ, eP=p−Shp+Shp−ph=θP+δP.

Let in (13). Then, we obtain:

 ρ2Δt(∥δn+1U∥2L2(Ω