Comput Geosci

DOI 10.1007/s10596-014-9418-y

ORIGINAL PAPER

On classical iterative subdomain methods for the Stokes–Darcy problem

Alfonso Caiazzo · Volker John · Ulrich Wilbrandt

Received: 18 December 2013 / Accepted: 20 March 2014 © Springer International Publishing Switzerland 2014

Abstract Within classical iterative subdomain methods, the problems in the subdomains are solved alternately by only using data on the interface provided from the other subdomains. Methods of this type for the Stokes–Darcy problem that use Robin boundary conditions on the interface are reviewed. Their common underlying structure and their main differences are identified. In particular, it is clarified that there are different updating strategies for the interface conditions. For small values of fluid viscosity and hydraulic permeability, which are relevant in applications from geosciences, it is shown in numerical studies that only one of these updating strategies leads to an efficient numerical method, if it is used with appropriate parameters in the

Robin conditions.

Keywords Stokes–Darcy problem · Classical iterative subdomain methods · Robin boundary conditions · Finite element methods · Continuous and discontinuous updating strategy · Applications from geosciences 1 Introduction

Consider a bounded domain Ω ⊂ Rd , d ∈ {2, 3}, and a decomposition Ω = Ωf ∪ ΓI ∪ Ωp into two disjoint subdomains Ωf and Ωp, denoting a free flow domain and a porous

A. Caiazzo · V. John () · U. Wilbrandt

Weierstrass Institute for Applied Analysis and Stochastics (WIAS), Mohrenstr. 39, 10117 Berlin, Germany e-mail: john@wias-berlin.de

V. John

Department of Mathematics and Computer Science, Free

University of Berlin, Arnimallee 6, 14195 Berlin, Germany medium, respectively, and possessing a common interface ΓI, i.e., Ωf ∩Ωp = ∅ and Ωf∩Ωp = ΓI. Assuming a moderate flow velocity in the free flow domain, the fluid dynamics in Ωf can be modeled with the incompressible Stokes equations for the velocity uf : Ωf → Rd [m/s] and the pressure

Pf : Ωf → R [Pa] ∇ · T ( uf,

Pf ρ ) = ff in Ωf, (1) ∇ · uf = 0 in Ωf. (2)

In Eq. (1), ρ [kg/m3] represents the fluid density, and the fluid stress tensor T(uf, pf) [m2/s2] is defined by

T(uf, pf) : Ωf → Rd×d, T = 2νD(uf) − pfI, where D(uf) = (∇u + ∇uT )/2 : Ωf → Rd×d [1/s] is the velocity deformation tensor, pf = Pf/ρ [m2/s2], and ν : Ωf → R [m2/s] denotes the kinematic viscosity of the fluid.

Outer forces acting on the free flow are modeled by ff : Ωf → Rd [m/s2].

The dynamics in the porous medium is described by the

Darcy law:

K∇ϕp + up = 0 in Ωp, (3) ∇ · up = fp in Ωp, (4) in terms of ϕp : Ωp → R [m], called Darcy pressure or piezometric head, and of the Darcy velocity up : Ωp → Rd [m/s]. In Eq. (3), K : Ωp → Rd×d [m/s] is the hydraulic conductivity tensor. Generally, it is assumed that K = KT ,

K > 0. Here, only the case K = KI, K > 0, will be

Comput Geosci considered. In Eq. (4), the function fp : Ωp → R [1/s] describes sinks and sources. System (3)–(4) is called the mixed form of the Darcy problem. An alternative formulation is obtained by taking the divergence of Eq. (3) and using (4), −∇ · (K∇ϕp) = fp in Ωp, (5) which is called the primal form of the Darcy problem.

The mixed form is often more important for applications, as it recovers the Darcy velocity up directly. However, we will focus on the primal formulation because of its simplicity. For the studied numerical methods, there are still open questions for this formulation.

System (1), (2), and (5) must be completed with appropriate boundary conditions and with proper interface conditions at the Stokes–Darcy interface ΓI. Denote by nf and np the unit outward normal vectors on ∂Ωf and ∂Ωp, respectively, and by τ i , i = 1, . . . , d − 1, pairwise orthogonal unit tangential vectors on the interface ΓI. Note that nf = −np on ΓI. Two standard coupling conditions on ΓI model the conservation of mass and the balance of normal stresses uf · nf = −K∇ϕp · nf on ΓI, (6) −nf · T(uf, pf) · nf = gϕp on ΓI, (7) where g [m/s2] is the gravitational acceleration. A classical third coupling condition is the so-called Beavers–Joseph–

Saffman condition [19, 25] (also Beavers–Joseph–Saffman–

Jones condition), which is based on experimental findings relating the tangential velocity along the interface to the fluid stresses: uf ·τ i+ατ i ·T(uf, pf)·nf = 0 on ΓI, i = 1, . . . , d−1, (8) where i = 1, . . . , d − 1, α = α0√(τ i · K · τ i )/(νg) [s/m], and α0 > 0 is a nondimensional parameter depending on properties of the porous medium. A simplification is based on the observation that both terms on the left-hand side of

Eq. (8) are small, leading to uf · τ i = 0 on ΓI, i = 1, . . . , d − 1. (9)

The coupled problem defined by Eqs. (1), (2), and (5) with the conditions (6), (7), and (8) (or (9)) has been studied extensively in the literature both theoretically; see, among others, [1, 6, 17, 21], and numerically, e.g., in [5, 9, 12–14, 16, 20, 24, 28].

Numerical methods for solving the coupled problem might be divided into two main classes. Direct, monolithic, or single-domain methods aim at solving the coupled system in a single step. Proposals of finite elementbased approaches in this class can be found, e.g., in [2–4, 16, 22, 26, 27]. As an alternative, decoupled, domaindecomposition, or multidomain approaches solve the coupled problem with a subdomain iterative procedure based, at each iteration, on the solution of the Stokes and

Darcy problem separately. On the one hand, an iterative method requires multiple solutions of the subproblems. On the other hand, these techniques allow to use specialized solvers for Stokes and Darcy problems and to tailor the algorithms according to their mathematical and physical properties, which might result in efficient procedures. Furthermore, iterative approaches are generally preferred if efficient solvers for the subproblems are available.