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.

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 [15]. 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 the -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 velocity norm on the Darcy subdomain.

2. Formulations of the Stokes-Darcy Coupled Equations

2.1. Model Equations

Let , or , be a bounded region with Lipschitz boundary . and are, respectively, the fluid and porous media subdomains of such that . The subdomains have a common interface . We denote by the fluid velocity and by the fluid pressure, where , , The flow in the domain is assumed to be of Stokes type and governed by the equations with symmetric strain tensor , external force , and constant viscosity In the porous region the filtration of an incompressible flow through porous media is described by Darcy equations where the permeability is a positive definite symmetric tensor and denotes an external Darcy force.

2.2. Boundary Conditions

On , we prescribe homogeneous Dirichlet conditions for the velocity . The boundary of is split into three parts We prescribe zero flux on and a homogeneous Dirichlet condition for the pressure on where denotes the outer normal vector on the boundary pointing from into . This boundary condition ensures a zero mass flux.

2.3. The Beavers-Joseph-Saffman Condition

The flows in and are coupled across the interface . Conditions describing the interaction of the flows are as follows [6, 8]:(i)The continuity of the normal velocity:(ii)The balance of normal forces:(iii)The Beavers-Joseph-Saffman condition written in terms of the strain tensor: where and is a dimensionless parameter to be determined experimentally, this condition relating the tangential slip velocity to the normal derivative of the tangential velocity component in the Stokes region

3. Variational Formulation

As variational formulation we consider the so-called - formulation used by Karper et al. [4] and recently by [9, 10]. We denote where and denote the usual Sobolev spaces.

Next, we define the spaces Then, multiplying the Stokes equations (1) by the test functions , , respectively, and integrating by part on the domain , we obtain Using the decomposition , the fluid normal stress condition (6), and the BJS interface condition (7) in (10), we obtain the weak formulation of the Stokes equations: find , such that ,

Similarly, taking and testing the Darcy equations (2) by , , respectively, together with the weighted grad-div term we obtain the weak formulation of Darcy equations: find , such that Summing up (11) and (12) the weak form of the coupled problem is given by the following: find , , , and such that

To analyze the weak formulation of the coupled problem we introduce the following spacesThe velocity and pressure spaces and are equipped with the natural norms

Further, due to the positive definiteness of with respect to the norm , there exist positive real numbers and such that

Next, we define the bilinear forms for , in and , in on the two parts of the domain by Hence, the bilinear form for the coupled problem is the sum of , , and terms to enforce the continuity of the normal part of the velocities across the interface. Assuming, for simplicity, that and are extended by zero to the whole domain, the variational formulation of the coupled Stokes-Darcy system in compact form reads as follows: find solution ofwithIt can easily be shown that a sufficiently regular solution of (19) such that , , is also a classical solution of (1) and (2). We note that there is an alternative variational formulation to the one given here called -formulation. The latter uses the term instead of [4].

The existence and uniqueness of the solution of problem (19) follows from Brezzi’s conditions for saddle point problems [11]; namely, with positive constants and [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 there is such that on , satisfying with positive constants , , and .

Proof. Let . Then, due to Stokes inf-sup condition there exists with on and on such that For the Darcy equation, due to the condition on , there exists with on and , such that Define and then where denote the Poincaré constant.
Then, using Young’s inequality we obtainChoosing , , positive constants such that we obtain the required result where In addition, we also havewhere .

4. Finite Element Discretization

Let be a shape-regular partition of quadrilaterals for or hexahedra for [12, 13]. The diameter of element will be denoted by and the global mesh size is defined by . Let be the reference element, the mapping from to element , and the space of all polynomials on with maximal degree in each coordinate. We assume that the mesh is obtained from a coarser mesh by global refinement. Hence, consists of patches of elements of We define the finite element space For the discrete spaces and we use the equal-order finite element functions that are continuous in and and piecewise polynomials of degree . We define the Scott-Zhang interpolation operator which preserves the boundary condition [13], as with stability and interpolation properties, respectively, aswhere , are positive constants.

We will also use the inverse inequality

Similarly, for vector functions we define the interpolation operatorwith 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 form by adding a stabilization term ; that is,such that the stabilized discrete problem reads 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 interpolant , onto a coarser mesh is used. Defining the subgrid stabilization term readswhere is patchwise constant.

A more attractive method from the computational point is obtained using only the fine mesh with smaller stencil. Defining the subgrid stabilization term reads

Next, we prove the stability of the discrete coupled Stokes-Darcy problem with respect to the norm

5. Stability

Theorem 2. Let be a quasi-regular partition [13]. Then, the following discrete inf-sup condition holds for some positive constant independent of the mesh size .

Proof. First, let , and then the diagonal testing combined with Korn’s inequality and the positivity of give In addition, let be as in Lemma 1, corresponding to , and set . Then,Next, we estimate and as follows:where the first two terms are bounded using Cauchy inequality together with the interpolation, stability, and inverse inequalitiesThe boundary term is bounded using the trace theorem and the - stability byHence, by Young inequality with we obtain where .
For the Darcy bilinear form we haveThen, by Young inequality and (52) we obtain Scaling we obtain Choosing we obtain which implies the required result with .

6. Error Analysis

Theorem 3. Assume that the solution of the Stokes-Darcy problem (19) is such that , , and is the solution of the stabilized problem (41). Then, the following error estimate holds with constants independent of :

Proof. Using the stability estimate of Theorem 3, there exists , with satisfying Then, by Galerkin orthogonality property, the first term of (59) is bounded by Hence, the approximation properties of and imply To estimate the second term of (59) we consider separately each individual term of the bilinear form .
Next, Cauchy schwarz and Poincaré inequality for the boundary terms imply Thus, Squaring the norm and applying Young inequality we obtain Next, we estimate the interpolation error by Adding the interpolation error (64) to the projection error (65) we obtain the required result

Remark 4. We note that the analysis above holds true for the triangular subgrid interpolation .

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 are and , respectively. So, we expect the -asymptotic rates to be and .

7. Numerical Results

As a test model problem we take and split it into and . The interface boundary is . We take , , , and and the right hand sides , such that the velocity and pressure solution in the two subdomains are given by 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: where , and .

The problem is solved using a velocity-pressure approximation with a two-level subgrid stabilization on a uniform mesh with . Rates of convergence for the velocity and pressure errors for , and are displayed in Tables 1 and 2.

In Table 1, we see clearly that the velocity field in the Stokes subdomain is of second-order accuracy with respect to the -norm and first-order accuracy with respect to -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 -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].

Competing Interests

The author declares that they have no competing interests.

Acknowledgments

The author acknowledges the financial support of the Sultan Qaboos University, under Contract IG/SCI/DOMS/14/07.