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  and a mortar multiscale finite element method presented in . 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 , 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 . 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 , 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 . Advection–diffusion on evolving diffuse interfaces has been studied in [33, 32]. We also mention the work in  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 . 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:
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:
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:
where is the unit normal to the fluid domain boundary.
The Beavers-Joseph-Saffman condition:
where , , is the orthonormal basis for the tangential space of and is the Beavers-Joseph-Saffman-Jones coefficient.
The balance of pressure:
Boundary conditions: We split the boundaries as: , , and prescribe the following boundary conditions:
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:
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 :
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  (see Section 3, conditions (S1)-(S3)). An example of such function is , where is a signed distance function, and
Then, we can write:
where is a Dirac distribution at the interface . Using the approximations above, and , we can approximate the interface integrals as:
where the tangent vectorsare obtained directly from the phase-field function using the algorithm described in [58, 31]. We introduce the following regularization of
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 , we define the weighted space associated to the measure by:
Associated with weighted spaces, we define the weighted Sobolev spaces:
with the corresponding norm
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:
For , we have , with equivalent norms given by:
Similarly, , with equivalent norms given by:
Notice that, due to the regularization (4), the phase field function is bounded as . Therefore, we have:
Because in the phase field formulation all unknowns are defined on the whole domain, , we define the following Hilbert spaces for :
Lemma 2 (Trace inequality)
Let be sufficiently small. Then, there exists a constant such that for and for , we have:
Lemma 3 (Poincaré inequality)
We define as the closure of in and set On these spaces, the following Poincaré inequality holds:
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 :
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.
There exists a unique weak solution in sense of Definition 2 to the Stokes-Darcy phase field formulation problem.
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
Notice that the solution is an element of . Now, we define a bilinear form on as:
and a linear functional as:
First we prove the following lemma.
The bilinear form is continuous and coercive on .
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:
Similarly, we have
Other terms can be bounded in the standard way. Therefore, we obtain the following estimate:
To show coercivity, we take to obtain:
Using the standard Korn and Poincare inequalities on together with Lemma 1, we have
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 :
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 . 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:
The discrete problem is given by: Find and such that for every and , we have:
3 Error analysis
We assume that for the chosen finite element spaces and , there exists a projection operator satisfying:
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  for and in  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
and let be a projection operator onto such that
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:
The main result of this section is stated in the following theorem.