Analysis of Subgrid Stabilization Method for Stokes-Darcy Problems

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.


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 [1][2][3][4][5]. 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 2 -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 1 (div) velocity norm on the Darcy subdomain.

Model Equations.
Let Ω ⊂ , = 2 or 3, be a bounded region with Lipschitz boundary Ω. Ω and Ω are, respectively, the fluid and porous media subdomains of Ω such that Ω ∩ Ω = 0. The subdomains have a common interface Γ = Ω ∩ Ω . We denote by k = (k , k ) the fluid velocity and by = ( , ) the fluid pressure, where k = k| Ω , = | Ω , = , . The flow in the domain Ω is assumed to be of Stokes type and governed by the equations −2] div ( (k )) + ∇ = f, in Ω div k = 0, in Ω (1) with symmetric strain tensor (k ) = (1/2)(∇k + ∇k ), external force f, and constant viscosity ] > 0. In the porous region Ω the filtration of an incompressible flow through porous media is described by Darcy equations 2 Advances in Numerical Analysis where the permeability = ( ) is a positive definite symmetric tensor and denotes an external Darcy force.

Boundary Conditions.
On Γ = Ω \ Γ, we prescribe homogeneous Dirichlet conditions for the velocity k .
where n denotes the outer normal vector on the boundary pointing from Ω into Ω . This boundary condition ensures a zero mass flux.

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: wherẽ= ] ⋅ and is a dimensionless parameter to be determined experimentally, this condition relating the tangential slip velocity k ⋅ to the normal derivative of the tangential velocity component in the Stokes region

3
To analyze the weak formulation of the coupled problem we introduce the following spaces The velocity and pressure spaces V and are equipped with the natural norms Further, due to the positive definiteness of with respect to the 2 (Ω ) norm ‖ ⋅ ‖ Ω , there exist positive real numbers 1 and 2 such that Next, we define the bilinear forms for k = (k , k ), w = (w , w ) in V and = ( , ), = ( , ) in on the two parts of the domain by Hence, the bilinear form for the coupled problem is the sum of A (k, ; w, ), A (k, ; w, ), and terms to enforce the continuity of the normal part of the velocities across the interface.
Assuming, for simplicity, that f 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 with It can easily be shown that a sufficiently regular solution is also a classical solution of (1) and (2). We note that there is an alternative variational formulation to the one given here called (div)-formulation. The latter uses the term −( , div w) The existence and uniqueness of the solution of problem (19) follows from Brezzi's conditions for saddle point problems [11]; namely, inf 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].
Proof. Let (k, ) ∈ . Then, due to Stokes inf-sup condition there exists w ∈ 1 (Ω ) with w = 0 on Γ and w ⋅ n = 0 on Γ such that For the Darcy equation, due to the condition = 0 on Γ ,2 , there exists w ∈ 1 (Ω ) with w ⋅ n = 0 on Γ ,2 and Γ, such that Define Advances in Numerical Analysis and then where denote the Poincaré constant. Then, using Young's inequality we obtain Choosing 1 , 2 , 3 positive constants such that we obtain the required result where In addition, we also have where 2 3 = max{ 2 , 2 ( 2 + 1)}.

Finite Element Discretization
Let T ℎ be a shape-regular partition of quadrilaterals for = 2 or hexahedra for = 3 [12,13]. The diameter of element ∈ T ℎ will be denoted by ℎ and the global mesh size is defined by ℎ fl max{ℎ : in T ℎ }. Let̂fl (−1; 1) be the reference element, the mapping from̂to element , and (̂) the space of all polynomials on̂with maximal degree ≥ 0 in each coordinate. We assume that the mesh T ℎ is obtained from a coarser mesh T 2ℎ by global refinement. Hence, T 2ℎ consists of patches of elements of T ℎ . We define the finite element space For the discrete spaces V ℎ and ℎ we use the equal-order finite element functions that are continuous in Ω and Ω and piecewise polynomials of degree ≥ 1.
We define the Scott-Zhang interpolation operator which preserves the boundary condition [13], as ℎ : 1 (Ω) → ℎ with stability and interpolation properties, respectively, as where , are positive constants. We will also use the inverse inequality Similarly, for vector functions we define the interpolation operator with 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 Advances in Numerical Analysis 5 on to the inf-sup condition. One possibility to circumvent this condition is to work with a modified bilinear form A ℎ (⋅; ⋅) by adding a stabilization term S ℎ (⋅; ⋅); 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 2ℎ , onto a coarser mesh T 2ℎ is used. Defining 2ℎ = − 2ℎ the subgrid stabilization term reads where 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

Stability
Theorem 2. Let T ℎ 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 (k ℎ , ℎ ) ∈ V ℎ × ℎ , and then the diagonal testing combined with Korn's inequality and the positivity of −1 give In addition, let w be as in Lemma 1, corresponding to (k ℎ , ℎ ) ∈ V ℎ × ℎ , and set z = j ℎ w − w. Then, Next, we estimate A (k ℎ , ℎ ; z, 0) and A (k ℎ , ℎ ; z, 0) as follows: where the first two terms are bounded using Cauchy inequality together with the interpolation, stability, and inverse inequalities The boundary term is bounded using the trace theorem and the 1 -stability by Hence, by Young inequality with we obtain

Theorem 3. Assume that the solution (v, ) of the Stokes-Darcy problem (19) is such that
is the solution of the stabilized problem (41). Then, the following error estimate holds with constants 1 , 2 , . . . , 7 independent of ℎ: Proof. Using the stability estimate of Theorem 3, there exists Then, by Galerkin orthogonality property, the first term of (59) is bounded by Advances in Numerical Analysis 7 Hence, the approximation properties of 2ℎ and ℎ imply To estimate the second term of (59) we consider separately each individual term of the bilinear form (1/̃)A ℎ (j ℎ k − k, ℎ − ; w ℎ , ℎ ). Next, Cauchy schwarz and Poincaré inequality for the boundary terms imply Thus, ) ℎ ‖k‖ +1,Ω + (̃3ℎ +̃4 ) ℎ ‖k‖ +1,Ω + (̃5 +̃6 1/2 ℎ 1/2 +̃7ℎ) ℎ +1,Ω 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 2 -asymptotic rates to be (ℎ +1 ) and (ℎ +1 ).
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 2 -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].