A-posteriori error analysis to the exterior Stokes problem

Applied Numerical Mathematics 63 (2013) 25–44

Contents lists available at SciVerse ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

A-posteriori error analysis to the exterior Stokes problem

✩

Mauricio A. Barrientos a,∗ , Matthias Maischak b a b

Instituto de Matemática, Pontificia Universidad Católica de Valparaíso, Blanco Viel 596, Cerro Barón, Valparaíso, Chile BICOM, Brunel University, Uxbridge, UB8 3PH, UK

a r t i c l e

i n f o

Article history: Received 17 August 2010 Received in revised form 9 May 2011 Accepted 14 September 2012 Available online 20 September 2012 Keywords: Exterior Stokes problem Dual-mixed finite element A-posteriori error estimates

a b s t r a c t We analyze the coupling of the dual-mixed finite element method with the boundary integral equation method. The result is a new mixed scheme for the exterior Stokes problem. The approach is based on the introduction of both the flux and the strain tensor of the velocity as further unknowns, which yields a two-fold saddle point problem as the resulting variational formulation. We show existence and uniqueness of the solution for the continuous and discrete formulations and provide the associated error analysis. In particular, the corresponding Galerkin scheme is defined by using piecewise constant functions and Raviart–Thomas spaces of lowest order. Most of our analysis makes use of an extension of the classical Babuška–Brezzi theory to a class of saddle-point problems. Also, we develop a-posteriori error estimates (of Bank–Weiser type) and propose a reliable adaptive algorithm to compute the finite elements solutions. Finally, several numerical results are given. © 2012 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Flow problems appear in many areas of science and engineering. One of the most important models is the Navier–Stokes problem. The non-linearity of the Navier–Stokes equations leads to difficulties in their numerical solution. As a consequence various simplifications of the Navier–Stokes equations have been investigated, cf. [25]. The Stokes equations are known as a model for a stationary incompressible viscous fluid with a small Reynolds number, which is the case for either high viscosity or very low velocity, the so-called creeping flows. For finite domains the Stokes equations are a good model in two and three-dimensions, whereas for infinite domains in 2d the Oseen’s equations are a more appropriate model. Due to technical reasons we restrict ourself to the Stokes equations in 2d. Our approach is also valid in 3d, and requires only a higher implementational effort. Also the Stokes equations (or Stokes-like equations) appear as linearizations of the Navier–Stokes equations, when Newton’s method is applied to solve the original non-linear problem. Classical analysis and implementation of the Stokes equations on finite domains are usually based on the so-called mixed finite element methods, cf. [9]. The primal-mixed finite element methods are usually based on a pressure-velocity variational formulation, where the discrete pressure field and the discrete velocity filed have to satisfy a compatibility condition, the so-called inf–sup or Babuška–Brezzi condition. E.g. in the simplest case we need P2-elements for the velocity and P1-elements for the pressure. The fast solution of pressure-velocity formulations is an active field of research, see e.g. the fictitious domain approach [2],

✩

*

This research was partially supported by Fondecyt-Chile through the research project No. 1070952. Corresponding author. E-mail addresses: [email protected] (M.A. Barrientos), [email protected] (M. Maischak).

0168-9274/$36.00 © 2012 IMACS. Published by Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.apnum.2012.09.003

26

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

domain decomposition [12,13], iterative solution with Uzawa type schemes and preconditioning [30,36] and boundary condition splitting [33] and the references therein. In the following we will concentrate on the dual-mixed approach, which allows us to compute stresses with higher accuracy than in the primal-mixed pressure-velocity approach. In case of the Stokes (and Stokes-like) problem we introduce the strain tensor and the rotation of the flux as auxiliary unknowns. The corresponding variational formulation leads to a two-fold saddle point operator equation. The abstract theory developed in [16], a slight generalization of the Babuška–Brezzi theory, is applied to prove the well-posedness of the continuous and discrete schemes. It is shown that the stability of the Galerkin scheme only requires low-order finite element subspaces, i.e. Raviart–Thomas spaces of order zero to approximate the strain and piecewise constant functions to approximate the other unknowns in the FEM part. More complicated elements as the PEERS space or the MINI element or Crouzeix–Raviart element are not necessary to satisfy the inf–sup condition. There exists a series of papers which analyze this formulation, the most precise analysis can be found in [15,14,17,18,11,22, 23]. Using the two-fold saddle point approach also allows to deal with non-linear Stokes equations, i.e. quasi-Newtonian flow, see e.g. [27]. On the other hand, the use of boundary integral equations for solving boundary value problems of Stokes-like problems has a long history, an initial work applying a BEM–FEM procedure for the exterior Navier–Stokes problem was given by [37]. This procedure starts from the so-called one boundary integral approach introduced by [29] for the Laplace equation. Usually, the discrete problem is posed on a polygonal approximation of an auxiliary boundary. This strategy has the serious drawback that the approximation of the nearly singular boundary integrals by simple quadrature is rendered seriously (see [29,37,38]). A more efficient method has been recently proposed in [32]. They discretize the integral operators on the original boundary (the regular auxiliary interface) and give a complete error analysis for both, the continuous and the fully discrete Galerkin method. Another approach to the external Stokes problem is given by using the finite element method with proper artificial boundary conditions by means of the Fourier series, cf. [26,4,5]. Moreover, [31] propose a new method to solve the twodimensional exterior Stokes problem numerically, using a non-symmetric primal-mixed fem–bem coupling. In [34] the corresponding symmetric primal-mixed fem–bem coupling method is numerically analyzed. In the following we will derive first the symmetric dual-mixed fem–bem coupling method for the exterior Stokes problem with an a-priori estimate. Second, as the main part of this work, we will develop an a-posteriori error estimator for the exterior Stokes problem, and will prove its reliability and quasi-efficiency. We prove in [6] that one can combine the classical Bank–Weiser method from [3] with the ideas from [10], to derive an a-posteriori error estimate for the two-fold saddle point formulation of a linear–non-linear transmission problem in plane elastostatics. We prove here that it can also successfully applied to Stokes problems in unbounded domains. All of that, the dual-mixed variational formulation (using the integral equations of the first and second kind in the unbound domain) and a-posteriori error estimates, up to the authors knowledge, is the first time applied to the exterior Stokes problem. At first glance, dual-mixed fem–bem coupling schemes can be seen as cumbersome, due to the large number of additional unknowns introduced in the weak formulation. But this approach offers certain advantages. E.g. variables of special interest, like the stress can be approximated more accurately. The system is symmetric, therefore an iterative solver like MINRES can be easily applied. Only standard spaces of lowest order are used, therefore the implementation is straight forward as well as preconditioning. The approach also allows easy extension to a local non-linear flow-model, like quasiNewtonian flow. Here we extend the original exterior Stokes problem in two ways. First we allow non-homogenous boundary conditions at the interior boundary. Second, we allow an un-physical behavior at infinity to simplify the incorporation of the integral operators. The physical behavior at infinity is restored by choosing the parameter Σ = 0. The rest of the paper is organized as follows. The model problem is described in Section 2. The dual-mixed variational formulation is shown in Section 3, its unique solvability is proved too. Section 4 deals with the Galerkin approximation of the problem. Piecewise constant and linear functions, and Raviart–Thomas spaces of order zero are used. Also, we demonstrate that it leads to a stable and convergent finite element scheme. In Section 5, an a-posteriori error analysis based on the Bank–Weiser method is developed and provide reliable and quasi-efficient a-posteriori error estimates. Several numerical results are reported in Section 6, which show a behavior according to the estimates obtained before. Finally we give some notation and definition of the spaces we will use: ·,· stands for the duality pairing between [ H −1/2 (Γ )]2 and [ H 1/2 (Γ )]2 with respect to the [ L 2 (Γ )]2 -inner product. Also, we introduce the following spaces:



2×2

L 2 (Ω)

    τ11 τ12 := τ = : τi j ∈ L 2 (Ω), i , j ∈ {1, 2} ,

H(div; Ω) :=

τ21 τ22





2×2

τ ∈ L 2 (Ω)



2

: div τ ∈ L 2 (Ω)

and

H0 (div, Ω) := with inner products



τ ∈ H(div, Ω): τ n, c  = 0 ∀c ∈ R2



M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

r, τ [ L 2 (Ω)]2×2 :=

r : τ dx,

r, τ H(div;Ω) :=

Ω

27

r : τ dx +

Ω

div r · div τ dx, Ω

and induced norms

 τ [ L 2 (Ω)]2×2 :=

2

1/2 τ

i , j =1

2 i j  L 2 (Ω)

,

1/2  . τ H(div,Ω) := τ 2[ L 2 (Ω)]2×2 +  div τ 2[ L 2 (Ω)]2 2. The model problem Let Ω0 be a simply connected and bounded domain in R2 with Lipschitz boundary Γ0 . Let f ∈ L 2 (R2 − Ω¯ 0 ) be given with bounded support and g ∈ [ H 1/2 (Γ0 )]2 . According to our hypothesis on f, we can choose a convex bounded annular domain Ω (from now on fixed) containing supp f and bounded by Γ0 and another Lipschitz-continuous closed curved Γ ¯ 0. whose interior contains supp f ∪ Ω Then we consider the steady-state exterior Stokes problem

−ν u + ∇ p = f in R2 − Ω¯ 0 , ∇ · u = 0 in R2 − Ω¯ 0 , u = g on Γ0 , u = Σ · log |x| + O (1)

as |x| → ∞,  −1  p = O |x| as |x| → ∞,

(2.1)

with Σ ∈ R2 be a given vector with



Σ=

σ · n ds, Γ

and σ the stress tensor described below. In what follows, R2×2 is the space of square matrices of order 2 with real entries, I  is the identity matrix of R2×2 , and 2 given τ := (τi j ), σ := (σi j ) ∈ R2×2 , we use the notations tr(τ ) := τ11 + τ22 and τ : σ = i , j =1 τi j σi j . The deviator of tensor

τ is denoted by dev(τ ) := τ − 12 tr(τ )I. We remark that tr(dev(τ )) = 0. Also, let us introduce the stress tensor σ ,

σ (u, p ) = − pI + 2ν e(u) with the strain tensor e(u) = 12 (∇ u + ∇ u T ). ¯ 0 into Ω and R2 − (Ω¯ 0 ∪ Ω), we can split (2.1) and write it, equivalently, as a Due to our split of the domain R2 − Ω Stokes problem in the bounded domain Ω :

−ν u + ∇ p = f in Ω, ∇ · u = 0 in Ω, u = g on Γ0 ,

(2.2)

coupled with the exterior and homogeneous Stokes problem:

−ν u + ∇ p = 0 in R2 − (Ω¯ ∪ Ω¯ 0 ), ∇ · u = 0 in R2 − (Ω¯ ∪ Ω¯ 0 ), u = Σ log |x| + O (1)

as |x| → ∞,  −1  p = O |x| as |x| → ∞,

(2.3)

by means of the following transmission conditions:

lim u =

x→x0 x∈Ω

lim

x→x0 x∈Ω

lim

x→x0

u,

¯ Ω¯ 0 ) x∈R −(Ω∪ 2

σ ·n=

lim

x→x0

¯ Ω¯ 0 ) x∈R −(Ω∪ 2

σ ·n

(2.4)

28

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

for almost all x0 ∈ Γ , where σ · n := σ (u, p ) · n denotes the stress vector at the boundary Γ , with n the unit outward normal to Ω in x0 . Also, with this notation it is possible to define the stress operators T (u) := σ (u, p ) · n and T (u) := σ (u, − p ) · n, which we will use to define the boundary integral operators in the next section. 3. The variational formulation The variational formulation arises from the coupling of boundary elements and mixed elements. Hence, in what follows ¯ 0) we use the weak formulation of (2.2)–(2.4) resulting from the application of the boundary integral method in R2 − (Ω ∪ Ω and the mixed finite element method in Ω . This formulation is obtained by incorporating the boundary integral equations satisfied by the solution of (2.3), into the dual-mixed variational formulation problem (2.2), by means of the transmission ¯ 0) conditions (2.4). These integral equations, which arise from the Green’s representation formula for u and σ in R2 − (Ω ∪ Ω and the jump conditions of the layer potentials, are given by



u= with

1 2



I + K u − V(σ n) + α

on Γ

(3.1)

α ∈ R2 , and   1 Wu + I + K (σ n) = 0 on Γ,

(3.2)

2

where V : [ H −1/2 (Γ )]2 → [ H 1/2 (Γ )]2 , K : [ H 1/2 (Γ )]2 → [ H −1/2 (Γ )]2 , K : [ H −1/2 (Γ )]2 → [ H 1/2 (Γ )]2 and W : [ H 1/2 (Γ )]2 → [ H −1/2 (Γ )]2 are the boundary integral operators of the simple, double, adjoint of the double and hypersingular layer potentials, respectively. More precisely, if

E (x, y ) :=

  (x − y )(x − y )T − log |x − y |I + 4πν |x − y |2 1

represents the fundamental velocity tensor, then we have





Vμ(x) :=

Kλ(x) :=

E (x, y )μ( y ) ds y ,

Γ \{x}





Γ \{x}



K μ(x) :=



 T x E (x, y ) μ( y ) ds y ,

T

T y E (x, y )



Wλ(x) :=

Γ \{x}

λ( y ) ds y ,

 

T

T x T y E (x, y )

Γ \{x}

λ( y ) ds y . 1 0

It is well known that (3.1) and (3.2) are not uniquely solvable [28], due to ker(V) = span{n} and ker(W) = span{ 0 ,  x2  −x } =: span{v1 , v2 , v3 } (rigid body motions). Therefore we have to replace (3.1) and (3.2) by (see [28, Eq. (2.3.36)]) 1



u=

1 2

1

,



V(σ n) + α I + K u −

on Γ

(3.3)

with the constraint



σ · n ds = Σ

(3.4)

Γ

(for given Σ ∈ R2 ) and

u+ W



1 2



I + K (σ n) = 0 on Γ,

(3.5)

with (see [28, Table 2.3.4])

 Vϕ := Vϕ +



 := Wu + ϕ n ds n, and Wu

Γ

3



uv j ds v j

on Γ.

(3.6)

j =1 Γ

As a first approach, we derive the weak formulation for the interior problem (2.2). To do that, we follow [17] regarding the value on the boundary. Thus, we use the strain tensor as an auxiliary unknown

t := e(u)

in Ω,

(3.7)

whence, the stress tensor has now the expression

σ = − pI + 2ν t in Ω.

(3.8)

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

29

In addition, since div(u) = tr(t) in Ω , we can write the incompressibility condition as

tr(t) = 0 in Ω.

(3.9)

According to the definition of e(u) we can rewrite (3.7) as

t = ∇u − γ , where

γ :=

(3.10)

1 (∇ u − ∇ u T ) 2





represents rotations and lies in the space

2×2

H 0 := δ ∈ L 2 (Ω)



: δ + δT = 0 =

Next, we multiply (3.10) by a test function



−



t : τ dx − Ω



u · div τ dx − Ω



0

δ

−δ 0



 : δ ∈ L 2 (Ω) .

(3.11)

τ ∈ H(div, Ω) and integrating by parts, we get

γ : τ dx + τ n, u∂Ω = 0.

(3.12)

Ω

Also, (3.8) and (3.9) yield, respectively,





t : s dx −

2ν Ω

and





Ω

2×2

p tr(s) dx = 0 ∀s ∈ L 2 (Ω)

σ : s dx −

(3.13)

,

Ω

q tr(t) dx = 0 ∀q ∈ L 2 (Ω),

−

(3.14)

Ω

and the equilibrium equation becomes



−



Ω



2

f · v dx ∀v ∈ L 2 (Ω) .

v · div σ dx =

(3.15)

Ω

Then, collecting (3.12), (3.13), (3.14), and (3.15), and requiring the symmetry of



σ weakly by

σ : δ dx = 0 for all δ ∈ H 0

(3.16)

Ω

we get the following problem: Find (t, p , σ , u, γ , ξ ) ∈ [ L 2 (Ω)]2×2 × L 2 (Ω) × H(div, Ω) × [ L 2 (Ω)]2 × H 0 × R such that





t : s dx −

2ν Ω

Ω

−

Ω





q tr(t) dx − Ω



p tr(s) dx = 0, Ω

t : τ dx −

−



σ : s dx −

Ω

Ω



v · div σ dx + η

u · div τ dx −

Ω

Ω



tr(σ ) dx =



γ : τ dx + ξ

tr(τ ) dx = −τ n, uΓ + τ n, gΓ0 , Ω

f · v dx, Ω

σ : δ dx = 0,

(3.17)

Ω

for all (s, q, τ , v, δ, η) ∈ [ L 2 (Ω)]2×2 × L 2 (Ω) × H(div, Ω) × [ L 2 (Ω)]2 × H 0 × R.  It is important to mention that in order to guarantee uniqueness we have required that Ω tr(σ ) dx = 0, which leads to the introduction of the new unknown ξ ∈ R and the test function η ∈ R. Now, we come back to Eqs. (3.3)–(3.6), which define the exterior problem on the boundary Γ , but before we couple it to the problem (3.13) we will give a more convenient way to represent these equations. Thus, introducing

 2 φ := u|Γ ∈ H 1/2 (Γ )

(3.18)

as further unknown, we can find an expression for φ from (3.5), as:

   −1 1 I + K (σ n), φ = −W 2

(3.19)

30

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

 : [ H 1/2 (Γ )]2 → [ H −1/2 (Γ )]2 being continuously invertible (the proof is similar to the Lamé case, see [28]). thanks to W Next, replacing (3.19) in (3.3), we get

    1 −1 1   φ = −V(σ n) − I+K W I + K (σ n) + α . 2

(3.20)

2

Whence, let R : [ H −1/2 (Γ )]2 → [ H 1/2 (Γ )]2 be the inverse Poincaré–Steklov operator given by

R (ϕ ) :=  V(ϕ ) +



1 2



 −1 I+K W



1 2



I + K (ϕ ).

(3.21)

Then

φ = − R (σ n) + α .

(3.22)

Consequently, if we replace (3.22) in the second equation of (3.17), we have: Find (t, p , σ , u, γ , ξ, α ) ∈ [ L 2 (Ω)]2×2 × L 2 (Ω) × H(div, Ω) × [ L 2 (Ω)]2 × H 0 × R × R2 such that





t : s dx −

2ν Ω

Ω

−



σ : s dx − Ω

t : τ dx −

Ω

p tr(s) dx = 0,



q tr(t) dx − Ω

u · div τ dx −

Ω



γ : τ dx + ξ Ω

tr(τ ) dx + τ n, α  =





τ n, R (σ n) + τ n, gΓ0 ,

Ω

σ n, β = β · Σ,



− v · div σ dx + η tr(σ ) dx = f · v dx,

Ω

Ω

Ω

σ : δ dx = 0,

(3.23)

Ω

for all (s, q, τ , v, δ, η, β) ∈ [ L 2 (Ω)]2×2 × L 2 (Ω) × H(div, Ω) × [ L 2 (Ω)]2 × H 0 × R × R2 . In order to give a dual-mixed weak formulation of the original problem (2.1), we write back our weak formulation using for R the boundary integral equations given by (3.3) and (3.5). Therefore, the weak formulation reads as: Find (t, φ, p , σ , u, γ , ξ, α ) ∈ [ L 2 (Ω)]2×2 × [ H 1/2 (Γ )]2 × L 2 (Ω) × H(div, Ω) × [ L 2 (Ω)]2 × H 0 × R × R2 such that



t : s dx −

2ν Ω



σ : s dx − Ω

Ω

Ω



u · div τ dx −

−

p tr(s) dx −



γ : τ dx + ξ Ω



t : τ dx −

Ω

q tr(t) dx + Ω

τ n,



1 2





V(σ n) I + K φ −

  1  tr(τ ) dx + W φ + I + K (σ n), μ + τ n, α  = τ n, gΓ0 , 



2

Ω

σ n, β = β · Σ,



− v · div σ dx + η tr(σ ) dx = f · v dx,

Ω

Ω

Ω

σ : δ dx = 0,

(3.24)

Ω

for all (s, μ, q, τ , v, δ, η, β) ∈ [ L 2 (Ω)]2×2 × [ H 1/2 (Γ )]2 × L 2 (Ω) × H(div, Ω) × [ L 2 (Ω)]2 × H 0 × R × R2 . We now show that (3.24) can be written in the form of a two-fold saddle point operator equation. To this purpose, we define the spaces



2×2

X 1 := L 2 (Ω)

 2 × H 1/2 (Γ ) ,

M 1 := L 2 (Ω) × H(div, Ω),



2

M := L 2 (Ω)

× H 0 × R × R2 ,

and define the operators A1 : X 1 → X 1 , B1 : X 1 → M 1 , B 1 : M 1 → X 1 , B : M 1 → M , B : M → M 1 , D : M 1 → M 1 , and the functional F ∈ M as follows

  A1 (r, ψ), (s, μ) := 2ν



Ω

 ψ, μ r : s dx + W

∀(r, ψ), (s, μ) ∈ X 1 ,

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

31

 



  1 ∀(r, ψ) ∈ X 1 , ∀(ρ , τ ) ∈ M 1 , B1 (r, ψ), (ρ , τ ) := − r : τ dx − ρ tr(r) dx + τ n, I+K ψ 



2

Ω

Ω

    B 1 (ρ , τ ), (r, ψ) := B1 (r, ψ), (ρ , τ ) ∀(ρ , τ ) ∈ M 1 , ∀(r, ψ) ∈ X 1 ,



  B(ρ , τ ), (v, δ, η, β) := − v · div τ dx − τ : δ dx + η tr(τ ) dx + τ n, β ∀(ρ , τ ) ∈ M 1 , ∀(v, δ, η, β) ∈ M , Ω

Ω

Ω

    B (v, δ, η, β), (ρ , τ ) := B(ρ , τ ), (v, δ, η, β) ∀(v, δ, η, β) ∈ M , ∀(ρ , τ ) ∈ M 1 ,     D(ρ , τ ), (λ, κ ) := κ n,  V(τ n) ∀(ρ , τ ), (λ, κ ) ∈ M 1 , G (ρ , τ ) = τ n, gΓ0 ∀(ρ , τ ) ∈ M 1 and

F (v, η, δ, β) := f · v dx + β · Σ ∀(v, η, δ, β) ∈ M . Ω

Hereafter, [·,·] stands for the duality pairing indicated by the corresponding operators and functional. Also, from now on O denotes a generic null operator/functional. According to the definitions above, the variational formulation (3.24) can be set as the following two-fold saddle point operator equation: Find ((t, φ), ( p , σ ), (u, γ , ξ, α )) ∈ X 1 × M 1 × M such that

⎞ ⎛ ⎞ ⎞⎛ A1 B 1 O O (t, φ) ⎝ B1 −D B ⎠ ⎝ ( p , σ ) ⎠ = ⎝ G ⎠ . F O B O (u, γ , ξ, α ) ⎛

(3.25)

Theorem 3.1. There exists a unique solution ((t, φ), ( p , σ ), (u, γ , ξ, α )) ∈ X 1 × M 1 × M of (3.25). Proof. Given that A1 is clearly continuous and X 1 -elliptic and that D is positive semi-definite, then following Theorem 4.3 of [19], we have to prove that B and B1 satisfy the corresponding inf–sup conditions. Thus, we want to prove that:

  [B(q, τ ), (v, δ, η, β)]  γ1 (v, δ, η, β) M (q, τ ) M 1 (q,τ )∈ M 1 \{0} sup

∀(v, δ, η, β) ∈ M

(3.26)

and

  [B1 (r, ψ), (q, τ )]  γ2 (q, τ ) M 1 (r, ψ) X 1 (r,ψ)∈ X 1 \{0}

1, ∀(q, τ ) ∈ M

sup

 1 := L 2 (Ω) × Q, and Q is defined as follows: where M

 

Q := τ ∈ H(div, Ω): div τ = 0 and τ = τ T in Ω, tr(τ ) dx = 0 and τ n = 0 , Ω

in others words, Q is just the kernel of the operator B . In fact, for (3.26) let (q, τ ) ∈ M 1 . Thus, we have

[B(q, τ ), (v, δ, η, β)] [B(0, τ ), (v, δ, η, β)]  sup (q, τ ) M 1 τ H(div,Ω) (q,τ )∈ M 1 \{0} τ ∈H(div,Ω)\{0}    − Ω v · div τ dx − Ω τ : δ dx + η Ω tr(τ ) dx + τ n, β = sup . τ H(div,Ω) τ ∈H(div,Ω)\{0} sup

Now, we observe that H(div, Ω) =  H(div, Ω) + RI, where

 

 H(div, Ω) := τ ∈ H(div, Ω): tr(τ ) dx = 0 . Ω

Indeed, τ = τ˜ + cI, with τ˜ ∈  H(div, Ω) and c =

  B(0, τ ), (v, δ, 0, β) = −



1 2|Ω|



v · div τ dx − Ω



Ω tr(τ ) dx ∈ R. Then, it is easy to see that

τ : δ dx + τ n, β Ω

    = B(0, τ˜ ), (v, δ, 0, β) = B(0, τ˜ ), (v, δ, η, β) .

(3.27)

32

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

Thus, it follows that

[B(q, τ ), (v, δ, η, β)] [B(0, τ˜ ), (v, δ, η, β)]  sup ( q , τ ) τ˜ H(div,Ω) M1 (q,τ )∈ M 1 \{0} τ˜ ∈ H(div,Ω)\{0} sup

[B(0, τ ), (v, δ, 0, β)] τ H(div,Ω) τ ∈H(div,Ω)\{0}   − Ω v · div τ dx − Ω τ : δ dx + τ n, β = sup . τ H(div,Ω) τ ∈H(div,Ω)\{0} =

sup

Next, given v ∈ [ L 2 (Ω)]2 , we let z ∈ [ H 1 (Ω)]2 be the weak solution of the mixed boundary value problem:

−z = v in Ω,

z = 0 on Γ0 ,

∂z = β on Γ, ∂n

(3.28)

for which one can easily show that z[ H 1 (Ω)]2  C 1 {v[ L 2 (Ω)]2 + |β|}. Then, we set in Ω , τ 0 n = β on Γ and τ 0 H(div,Ω)  C 2 {v[ L 2 (Ω)]2 + |β|}. It follows that

τ 0 = ∇ z and observe that − div τ 0 = v

[B(q, τ ), (v, δ, η, β)] [B(0, τ 0 ), (v, 0, 0, β)]  sup (q, τ ) M 1 τ 0 H0 (div,Ω) (q,τ )∈ M 1 \{0} τ 0 ∈H(div,Ω)\{0} sup

= where

v2[ L 2 (Ω)]2 + |Γ ||β|2 τ 0 H(div,Ω)

  γ1 v[ L 2 (Ω)]2 + |β| ,

(3.29)

γ1 depends on |Γ | and C 2 . Also, we know that (see Theorem 4.3 of [19] or Lemma 4.3 in [8]) there exists   [B(0, τ ), (v, δ, 0, 0)]  γ¯ (v, δ)[ L 2 (Ω)]2 × H . 0 τ H(div,Ω) τ ∈H(div,Ω)\{0} sup

Finally, noting that

1/2 [B(q, τ ), (v, δ, η, β)] [B(0, ηI), (v, δ, η, 0)]   = 2|Ω| |η|, ( q , τ )  η I  M1 H(div,Ω) (q,τ )∈ M 1 \{0} sup

(3.30)

the proof of (3.26) is completed.  1 and assume that q L 2 (Ω)  τ H(div,Ω) . Thus, Now, for the case of (3.27) let (q, τ ) ∈ M

sup

(r,ψ)∈ X 1 \{0}

[B1 (r, ψ), (q, τ )] [B1 (r, 0), (q, τ )]  sup (r, ψ) X 1 2 2 × 2 r∈[ L (Ω)] \{0} r[ L 2 (Ω)]2×2   − Ω r : τ dx − Ω q tr(r) dx = sup . r[ L 2 (Ω)]2×2 r∈[ L 2 (Ω)]2×2 \{0}

Then, we can choose r = −[τ −

1 2

tr(τ )I] (which has tr(r) = 0!), and following the proof of Lemma 3.3 of [24] we have,

[B1 (r, ψ), (q, τ )]  τ [ L 2 (Ω)]2×2 = τ H(div,Ω) . (r, ψ) X 1 (r,ψ)∈ X 1 \{0} sup

On the other hand, if q L 2 (Ω)  τ H(div,Ω) , then taking r = −qI + τ , we get again

[B1 (r, ψ), (q, τ )] [B1 (r, 0), (q, τ )]  sup (r, ψ) X 1 r[ L 2 (Ω)]2×2 (r,ψ)∈ X 1 \{0} r∈[ L 2 (Ω)]2×2 \{0} sup

 which completes the proof of (3.27).

2q2L 2 (Ω) − τ 2H(div,Ω)

qI + τ [ L 2 (Ω)]2×2 2



1 1+

√ qL 2 (Ω) , 2

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

33

4. A finite element solution

L

Hereafter, we assume that Γ and Γ0 are polygonal curves. Also, let us assume that the boundary Γ = j =1 Γ j consists of L separate, mutually non-intersecting compact boundary components Γ1 , . . . , Γ L . Then, we let{Th }h>0 be a regular family ¯ = { T : T ∈ Th }. Also, for each of triangulations of Ω , made up of triangles T of diameter h T , such that h := sup h T and Ω T ∈Th

T ∈ Th we let R T 0 ( T ) be the local Raviart–Thomas space of order zero, that is

R T := span

      1 0 x , , 1 , 0

1

x2

where (x1 , x2 ) ∈ R is a generic vector, [35]. In addition, given a non-negative integer k and S ⊂ R2 , we define the space of polynomials Pk (S ) on S of degree  k. Thus, we define the following finite elements subspaces T

2





2×2

X ht := sh ∈ L 2 (Ω) φ

X h :=







2×2

: sh | T ∈ P0 ( T )

2



2

∀ T ∈ Th ,



μh ∈ H 1/2 (Γ ) : μh |T ∈ P1 (Γ j ) , j = 1, . . . , L , φ

X 1,h := X ht × X h ,





p

M h := q ∈ L 2 (Ω): q| T ∈ P0 ( T ) ∀ T ∈ Th ,



M hσ :=





τ := (τi j ) ∈ H(div, Ω): (τi1 , τi2 )T T ∈ R T 0 ( T ) ∀i ∈ {1, 2}, ∀ T ∈ Th , p

M 1,h := M h × M hσ ,





2



2

M hu := vh ∈ L 2 (Ω) : vh | T ∈ P0 ( T )



H 0h := M h :=

M hu

0

δ

−δ 0



∀ T ∈ Th ,

  2×2 ∈ H 1 (Ω) : δ| T ∈ P 1 ( T ) ∀ T ∈ Th ,

× H 0h × R × R2 .

Consequently, the Galerkin scheme associated with (3.25) reads: Find ((th , φh ), ( p h , σ h ), (uh , γ h , ξh , αh )) ∈ X 1,h × M 1,h × M h such that

⎞ ⎛ ⎞ ⎞⎛ A1 B 1 O O (th , φh ) ⎠ = ⎝ G ⎠. ⎝ B1 −D B ⎠ ⎝ ( ph , σ h ) F O B O (uh , γ h , ξh , αh ) ⎛

(4.1)

Theorem 4.1. The discrete problem (4.1) is uniquely solvable and the following error estimate holds:

  (t, φ, p , σ , u, γ , ξ, α ) − (th , φh , ph , σ h , uh , γ , ξh , αh ) h H    C inf (t, φ, p , σ , u, γ , ξ, α ) − (sh , μh , qh , τ h , vh , δh , ηh , βh )H , h

(sh ,μh ,qh ,τ h ,vh ,δh ,ηh ,βh )∈Hh

where H := X 1 × M 1 × M, Hh := X 1,h × M 1,h × M h and C is a constant independent of h. Proof. As in Theorem 3.1, we note again that A1 is clearly continuous and X 1,h -elliptic and that D is positive semi-definite, then applying again the abstract theory of [19], we have to prove that B and B1 satisfy the corresponding discrete inf–sup conditions with constants independent of h. To do that, we proceed as Theorem 3.1 of [17]. Given (vh , δ h , ηh , βh ) ∈ M h , we have the discrete analogue of (3.30), that is

1/2 [B(qh , τ h ), (vh , δ h , ηh , βh )] [B(0, ηh I), (vh , δ h , ηh , 0)]   = 2|Ω| |ηh |. (qh , τ h ) M 1 ηh IH(div,Ω) (qh ,τ h )∈ M 1 \{0} sup

ˆ σ + RI, where M ˆ σ := M σ ∩  Now, we observe that M hσ = M H(div, Ω). Thus, it follows that h h h sup

(qh ,τ h )∈ M 1,h \{0}

[B(qh , τ h ), (vh , δh , ηh , βh )] [B(0, τ˜ h ), (vh , δh , ηh , βh )]  sup (qh , τ h ) M 1 τ˜ h H(div,Ω) σ ˆ τ˜ h ∈ M h \{0} [B(0, τ h ), (vh , δh , 0, βh )] τ h H(div,Ω) τ h ∈ M hσ \{0}   − Ω vh · div τ h dx − Ω δh : τ h dx + τ h n, βh  = sup . τ h H(div,Ω) τ h ∈ M σ \{0} =

sup

h

34

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

Then, given (vh , δ h ) ∈ M hu × H 0h , we know (see Lemma 4.4 in [1]) that there exists

τ h ∈ M hσ , τ h = 0, such that

    B(0, τ h ), (vh , δ h , 0, 0)  γ¯ τ h H(div,Ω) (vh , δ h )[ L 2 (Ω)]2 × H , 0h

with the constant γ¯ independent of h. On the other hand, we have

− [B(0, τ h ), (vh , 0, 0, βh )] sup  sup σ  τ  τ h,i ={0} h H(div,Ω) τ h ∈ M h \{0}



Ω v h,i · div τ h,i dx + τ h,i n, βh 

τ h,i H(div,Ω)

where v h,i and τ h,i are the i-th component and i-th row of the vector vh and tensor τ h , respectively. Then, by regularity reasons we can’t get directly results like (3.28) and (3.29), but following closely the proof given in Lemma 3.4 of [7], it can be proved that

sup τh

∈ M σ \{0} h

 [B(0, τ h ), (vh , 0, 0, βh )]  γ ∗  v h,i L 2 (Ω) + |β| , τ h H(div,Ω)

which proves the discrete inf–sup condition for B .  1,h := M p × Qh , where Qh is defined as: Further, let M h

 

Qh := τ h ∈ M hσ : div τ h = 0 in Ω, tr(τ h ) dx = 0 , Ω

the discrete kernel of B . Then, to prove the discrete inf–sup condition for B1 , it is necessary to distinguish two cases: when qh  L 2 (Ω)  τ h H(div,Ω) and qh  L 2 (Ω)  τ h H(div,Ω) . In both cases, the proof is quite similar to Theorem 3.1 above or Lemma 3.3 of [24]. Finally, following the proof of Theorems 6.1 and 5.1 of [19] the proof is completed. 2 5. A-posteriori error analysis In this section, we follow the approach from [6] (see also [24]) to derive a Bank–Weiser’s type a-posteriori estimate for the error

  (t, φ, p , σ , u, γ , ξ, α ) − (th , φh , ph , σ h , uh , γ , ξh , αh ) . h H For this purpose, let X := X 1 × M 1 and introduce the linear saddle point operator A :=

 A1

B1  :X B1 −D

      A(r, ψ, λ, κ ), (s, μ, ρ , τ ) := A1 (r, ψ), (s, μ) + B1 (r, ψ), (ρ , τ )     + B1 (r, ψ), (ρ , τ ) − D(ρ , τ ), (λ, κ ) ,

→ X , that is

(5.1)

for all (r, ψ, λ, κ ), (s, μ, ρ , τ ) ∈ X . In addition, let Aˆ : X → X be the linear and bounded operator induced by the inner product of X , that is

   ψ, μ + λ, ρ  L 2 (Ω) + κ , τ H(div,Ω) , Aˆ (r, ψ, λ, κ ), (s, μ, ρ , τ ) := r, s[ L 2 (Ω)]2×2 + W

(5.2)

for all (r, ψ, λ, κ ), (s, μ, ρ , τ ) ∈ X . ¯ p¯ , σ¯ ) ∈ X Thus, we can define the X -Ritz projection of the finite element error, with respect to Aˆ , as the unique (t¯, φ, such that

    ¯ p¯ , σ¯ ), (s, μ, q, τ ) := A(t − th , φ − φh , p − ph , σ − σ h ), (s, μ, q, τ ) Aˆ (t¯, φ,   + B(q, τ ), (u, γ , ξ, α ) − (uh , γ h , ξh , αh )

(5.3)

¯ p¯ , σ¯ ) is guaranteed by the fact that the right-hand side of (5.3) constitutes a for all (s, μ, q, τ ) ∈ X . The existence of (t¯, φ, linear and bounded functional in X . Next, given T ∈ Th we denote by ·,·[ L 2 ( T )]2×2 , ·,· L 2 ( T ) , and ·,·H(div, T ) the inner products of [ L 2 ( T )]2×2 , L 2 ( T ), and H(div, T ), respectively, and let (th, T , p h, T , σ h, T , uh, T , γ h, T ) be the restriction of (th , p h , σ h , uh , γ h ) to T .

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

35

5.1. Main results The following theorem provides an important upper bound for the Ritz projection of the Galerkin error with respect to the operator Aˆ .

˜ h be a function in [ H 1 (Ω)]2 such that ϕ˜ h (¯x) = g(¯x) for each vertex x¯ of Th lying on Γ0 . For each T ∈ Th define Theorem 5.1. Let ϕ tˆ T := σ h, T − (2ν th, T − p h, T I),

(5.4)

and let σˆ T ∈ H(div, T ) be the unique solution of the local problem

σˆ T , τ H(div,T ) = F h,T (τ ) ∀τ ∈ H(div, T ), where F h, T ∈ H(div, T ) is defined by



F h, T (τ ) :=





(th,T + γ h,T ) : τ dx + T

(5.5)

T



2



ϕ˜ h · τ nT ds +

tr(τ ) dx − T

  1 + τ n, ϕ˜ h +  V(σ h n) − I + K φh − αh 



uh, T · div τ dx − ξh

(g − ϕ˜ h ) · τ nT ds

∂ T ∩Γ0

∂T

(5.6)

,

∂ T ∩Γ

with n T being the unit outward normal to ∂ T . Then

    tr(th,T )22 + ¯ p¯ , σ¯ ), (t¯, φ, ¯ p¯ , σ¯ )  Aˆ (t¯, φ, tˆ T 2[ L 2 (T )]2×2 + σˆ T 2H(div,T ) L (T ) T ∈Th

T ∈Th

 2    −1     φh + 1 I + K (σ h n) +  W  W   2

T ∈Th

[ H −1/2 (Γ )]2

(5.7)

.

Proof. From the variational formulation (3.25) we deduce that



    A(t, φ, p , σ ), (s, μ, q, τ ) + B(q, τ ), (u, γ , ξ, α ) =

g · τ n ds,

(5.8)

Γ0

whence, (5.3) yields

  ¯ p¯ , σ¯ ), (s, μ, q, τ ) = Aˆ (t¯, φ,











gτ n ds − A(th , φh , p h , σ h ), (s, μ, q, τ ) − B (q, τ ), (uh , γ h , ξh , αh ) .

(5.9)

Γ0

Since Aˆ is symmetric and X -elliptic, we have that



1

 ¯ p¯ , σ¯ )2 = − (t¯, φ, X 2

min

(s,μ,q,τ )∈ X

    (s, μ, q, τ )2 − Aˆ (t¯, φ, ¯ p¯ , σ¯ ), (s, μ, q, τ ) , X

1 2

which, using (5.9), becomes

 1 ¯ p¯ , σ¯ )2 = − (t¯, φ, X 2

min

(s,μ,q,τ )∈ X

J(s, μ, q, τ )

where

J(s, μ, q, τ ) :=

1 2

 (s, μ, q, τ )2 − X

(5.10)

gτ n ds Γ0

    + A(th , φh , ph , σ h ), (s, μ, q, τ ) + B(q, τ ), (uh , γ h , ξh , αh ) .

(5.11)

˜ h belongs to [ H 1 (Ω)]2 , we can write Now, since the function ϕ

−





ϕ˜ h · τ nT ds +

T ∈Th ∂ T



ϕ˜ h · τ n ds + Γ

ϕ˜ h τ n ds = 0.

(5.12)

Γ0

Then, adding (5.12) to (5.11) and using the definitions of A and B , we obtain

J(s, μ, q, τ ) :=



T ∈Th

J1, T (s T ) + J2 (μ) +



T ∈Th

J3, T (q T ) +



T ∈Th

J4, T (τ T ),

(5.13)

36

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

where (s T , q T , τ T ) is the restriction of (s, q, τ ) to T ,

1

sT 2[ L 2 (T )]2×2 − tˆ T , sT [ L 2 (T )]2×2 ,     1 1   J2 (μ) := Wμ, μ + Wφh + I + K (σ h n), μ ,

J1, T (s T ) :=

2

J3, T (q T ) :=

(5.14)

2

(5.15)

2

1 2

  q T 2L 2 (T ) + tr(th,T ), q T L 2 (T ) ,

(5.16)

and

J4, T (τ T ) :=

1 2

τ T 2H(div,T ) − F h,T (τ T ).

(5.17)

Noting that tr(th, T ) ∈ L 2 ( T ) and hence it follows that

2 1 min J3, T (q T ) = − tr(th, T ) L 2 ( T ) . 2

q T ∈L2 (T )

2

The rest of the proof is straightforward to the proof of Theorem 3.1 of [6].

The following theorem establishes a first a-posteriori error estimate for our finite element solution. For this, we use the X -Ritz projection of the error and the corresponding upper bound in Theorem 5.1.

˜ h be a function in [ H 1 (Ω)]2 . For each T ∈ Th define Theorem 5.2. Let ϕ tˆ T := σ h, T − (2ν th, T − p h, T I), and let σˆ T ∈ H(div, T ) be the unique solution of the local problem (5.5). Then, there exists C > 0, independent of h, such that

  (t, φ, p , σ , u, γ , ξ, α ) − (th , φh , ph , σ h , uh , γ , ξh , αh )  C h H



1/2 θ 2T + R2Γ + R2Σ

,

(5.18)

T ∈Th

where

2  2  θ 2T := tˆ T 2[ L 2 (T )]2×2 + tr(th,T ) L 2 (T ) + σˆ T 2H(div,T ) + f + div σ h,T 2[ L 2 (T )]2 + σ h,T − σ hT,T [ L 2 (T )]2×2 ,       1   RΓ := Wφh + I + K (σ h n)  2

(5.19) (5.20)

[ H −1/2 (Γ )]2

and





RΣ := Σ − σ h n, 1.

(5.21)

Proof. Taking note that operator A1 is a uniformly bounded and X 1 -elliptic bilinear form, that B and B1 satisfy the continuous inf–sup conditions and that D is positive semi-definite, we conclude by applying Theorem 2.3 of [20] (see also Theorem 3.3 of [6]) that the linear operator obtained by adding the three equations of the left-hand side of (3.25) satisfies a global inf–sup condition, with a constant C 0 > 0 such that

  (t, φ, p , σ , u, γ , ξ, α ) − (th , φh , ph , σ h , uh , γ , ξh , αh ) h H      C0 sup A1 (t − th , φ − φh ), (s, μ) + B1 (t − th , φ − φh ), (q, τ ) (s,μ,q,τ ,v,δ,η,β)H 1

   + B1 (s, μ), ( p − ph , σ − σ h ) − D(q, τ ), ( p − ph , σ − σ h )     + B(q, τ ), (u − uh , γ − γ h , ξ − ξh , α − αh ) + B( p − ph , σ − σ h ), (v, δ, η, β) . 

(5.22)

Now, using the X -Ritz projection of the error (5.3), the third equation of (3.25) and the definitions of B and F , we obtain that

1 

 (t, φ, p , σ , u, ξ, α ) − (th , φh , ph , σ h , uh , ξh , αh ) H   



  ¯ p¯ , σ¯ ), (s, μ, q, τ ) + σ h : δ dx + (f + div σ h ) v dx + Σ − σ h n, β .  sup Aˆ (t¯, φ, |Γ | (s,μ,q,τ ,v,δ,η,β)H 1 C0

Ω

Ω

(5.23)

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

37

Next, by the boundedness and ellipticity of Aˆ and Theorem 5.1, we see that

       Aˆ (t¯, φ, ¯ p¯ , σ¯ ), (s, μ, q, τ )   Aˆ (t¯, φ, ¯ p¯ , σ¯ ) (s, μ, q, τ ) X X   12   2    2 2 2   ˆ tT [ L 2 (T )]2×2 + tr(th,T ) L 2 (T ) + σˆ T H(div,T ) + RΓ (s, μ, q, τ ) X , C T ∈Th

(5.24) for all (s, μ, q, τ ) ∈ X . Also, it is possible to see that



 

   12     1   T T  σ h : δ dx =  1  σ h − σ h : δ dx  σ h − σ h [ L 2 (T )]2×2 δ[ L 2 (Ω)]2×2 ,   2 2 Ω

(5.25)

T ∈Th

Ω



   12   2  (f + div σ h ) v dx  f + div σ h [ L 2 (T )]2 v[ L 2 (Ω)]2   Ω

and

(5.26)

T ∈Th

   Σ     RΣ |β|. σ n , β − h  Γ 

(5.27)

Therefore, replacing (5.24), (5.25), (5.26) and (5.27) into (5.23), we complete the proof.

2

5.2. Quasi-efficiency The next lemma give us an a-priori estimate for the solution of (5.5). Moreover, it will be used to study the efficiency of the estimate obtained in the previous theorem and to deduce a fully explicit a-posteriori error estimate for the finite element solution.

˜ h and σˆ T be as indicated in Theorem 5.1. Then, there exists C > 0, independent of h, ν and T , such that Lemma 5.3. Let ϕ

 ˆ σ T H(div,T )  C th,T + γ h,T − ∇ ϕ˜ h 2[ L 2 (T )]2×2 + uh,T − ϕ˜ h 2[ L 2 (T )]2 + g − ϕ˜ h 2[ H 1/2 (∂ T ∩Γ

0 )]

2      1   ˜ φ + ϕ + V ( σ n ) − + K − α I h h h   h 2

2

+ h2T |ξh |2

1/2 .

(5.28)

[ H 1/2 (∂ T ∩Γ )]2

In addition, for any z ∈ [ H 1 (Ω)]2 such that z = g on Γ0 , we get

 σˆ T H(div,T )  C th,T + γ h,T − ∇ z2[ L 2 (T )]2×2 + uh,T − z2[ L 2 (T )]2 + h2T |ξh |2 + z − ϕ˜ h 2[ H 1/2 (∂ T )]2 2      1   φ + z + V ( σ n ) − + K − α I h h h  2

1/2 ,

(5.29)

[ H 1/2 (∂ T ∩Γ )]2

where ∂ T is the part of ∂ T not contained in Γ ∪ Γ0 .

˜ h is a function in [ H 1 (Ω)]2 , a simple integration by parts gives Proof. Since ϕ







∇ ϕ˜ h : τ dx +

ϕ˜ h · τ nT ds = T

∂T

ϕ˜ h · div τ dx. T

Thus, replacing the above expression back into F h, T , applying Cauchy–Schwarz’s inequality, and remembering that σˆ T ∈ H(div, T ), as the Riesz representant of the functional F h, T ∈ H(div, T ) , is such that σˆ T H(div, T ) =  F h, T H(div, T ) , we obtain (5.28). Now, given z ∈ [ H 1 (Ω)]2 such that z = g on Γ0 , we can write



− ∂T



ϕ˜ h · τ nT ds + ∂ T ∩Γ0

(ϕ˜ h − g) · τ nT ds

38

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

=−



∂T

∂T



=−

(z − ϕ˜ h ) · τ nT ds +

z · τ n T ds +

∇ z : τ dx −

T

∂ T ∩Γ0



(z − ϕ˜ h ) · τ nT ds,

z · div τ dx + T

(ϕ˜ h − z) · τ nT ds (5.30)

∂T

which, replaced back into F h, T , gives (5.29) and concludes the proof.

2

The above lemma motivates the following result, which shows the quasi-efficiency of θ T given by Theorem 5.2, which ˜ h ) on the edges of Th in the interior of Ω . means that it is efficient up to a term depending on the traces (u − ϕ

˜ h ∈ [ H 1 (Ω)]2 . Then there exists C > 0, independent of h, such that for all T ∈ Th Lemma 5.4. Assume that u and ϕ



θ 2T

 C t − th,T 2[ L 2 (T )]2×2 + u − uh,T 2[ L 2 (T )]2 +  p − ph,T 2L 2 (T ) + σ − σ h,T 2H(div,T ) + h2T |ξ − ξh |2 + u − ϕ˜ h 2[ H 1/2 (∂ T )]2  2 + γ − γ h,T 2[ L 2 (T )]2×2 +  V(σ n − σ h n)[ H 1/2 (∂ T ∩Γ )]2  2    1  2  (φ − φ + + K ) +  α − α  I h  h [ H 1/2 (∂ T ∩Γ )]2 ,  2

and hence

T ∈Th

(5.31)

[ H 1/2 (∂ T ∩Γ )]2

  2 θ 2T + R2Γ + R2Σ  C (t, φ, p , σ , u, γ , ξ, α ) − (th , φh , ph , σ h , uh , γ h , ξh , αh )H +

e ∈ E h (Ω)

 u − ϕ˜ h 2[ H 1/2 (e)]2 ,

(5.32)

where E h (Ω) is the set of edges of Th in the interior of Ω . Proof. By taking s = μ = 0 in the first equation of (3.24), we have



−





t : τ dx − Ω

q tr(t) dx +



τ n,

1 2





 tr(τ ) dx = τ n, gΓ0 . I + K φ − V(σ n) + α − u · div τ dx + ξ

Ω

Ω

Ω

In addition, we easily get ξ = 0 and tr(t) = 0 in Ω . Then, since u ∈ [ H (Ω)] , we can integrate by parts in the above equation and it reduces to 1



2

      1  ∇ u − (t + γ ) : τ dx + τ n, I + K φ − V(σ n) − u + α + τ n, uΓ0 = τ n, gΓ0 . 2

Ω

It follows that ∇ u = t + γ in Ω , u = g on Γ0 and u = ( 12 I + K)φ − V(σ n) + α on Γ . Then, applying (5.29) with z = u, we deduce that

 σˆ T 2H(div,T )  2 th,T − t2[ L 2 (T )]2×2 + γ h,T − γ 2[ L 2 (T )]2×2

 2 + uh − u2[ L 2 (T )]2 + u − ϕ˜ h 2[ H 1/2 (∂ T )]2 + h2T |ξh − ξ |2 + V(σ h n − σ n)[ H 1/2 (∂ T ∩Γ )]2  2    1  2  (φ + + K − φ) +  α − α  I h h [ H 1/2 (∂ T ∩Γ )]2 .   2

(5.33)

[ H 1/2 (∂ T ∩Γ )]2

Also, we have that

  tr(th,T ) Now, taking

L2 (T )

  = tr(th,T − t) L 2 (T )  th,T − t[ L 2 (T )]2×2 .

τ = μ = 0 in the first equation of (3.25), we observe that

σ = 2ν t − pI in Ω,

(5.34)

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

39

and hence

  tˆ T [ L 2 (T )]2×2 = σ h,T − (2ν th,T − ph,T I)[ L 2 (T )]2×2    (σ h,T − σ ) + (2ν t − p I) − (2ν th,T − ph,T I)[ L 2 (T )]2×2  σ h,T − σ H(div,T ) + 2ν t − th,T [ L 2 (T )]2×2 + 2 p − ph,T L 2 (T ) . Next, it follows from the third equation of (3.25) that − div σ = f in Ω , Whence, we get

σ = σ T in Ω and

(5.35)



Ω tr(σ ) dx = 0.

  f + div σ h,T [ L 2 (T )]2 = div(σ − σ h,T )[ L 2 (T )]2 ,

(5.36)

  σ h,T − σ T  2 2×2  2σ h,T − σ  2 2×2 . [ L ( T )] h, T [ L ( T )]

(5.37)

and

On the other hand, also we find from (3.24) that

φ+ W



1 2



I + K (σ n) = 0,

 and K , implies that which, by the continuity of the boundary integral operators W





R2Γ  C φ − φh 2[ H 1/2 (Γ )]2 + σ − σ h 2H(div,Ω) .

(5.38)

Finally, the estimates (5.33), (5.34), (5.35), (5.36) and (5.37) imply the local quasi-efficiency of the error estimator θ 2T . The V and K, summing up over all the elements T ∈ Th and (5.38). 2 proof is completed by the continuity of the operators  At this point, we remark that the a-posteriori error estimate provided by Theorem 5.2 includes the non-local term RΓ , which is defined in terms of a negative order Sobolev norm. In addition, it is important to note that the local problem defining σˆ T H(div, T ) is set in a space of infinite dimension H(div, T ). This implies that (5.5) must be solved by using, for instance, the h, p or h − p version of the finite element method, which yields approximations of the local indicators θ T . Nevertheless, the main property of θ T , as proved in the above lemma, is that it constitutes a quasi-efficient and reliable a-posteriori error estimate. 5.3. Fully explicit a-posteriori error estimate Now, our purpose is to bound the contributions RΓ and σˆ T H(div, T ) by easily computable expressions. First, we start with the following estimation of RΓ : Lemma 5.5. There exists a constant C > 0, independent of h, such that

 2   L    1   RΓ  C log 1 − C h (Γ ) h j Wφh + I + K (σ h n)  

j =1

where h j = length of Γ j , and



C h (Γ ) := max

hi hj

2

(5.39)

,

[ L 2 (Γ j )]2



: | i − j | = 1, i , j ∈ {1, . . . , L } .

Proof. See Lemma 3.5 of [6] (see also Corollary 1 of [21]).

2

In consequence, we can stablish the a-posteriori error estimate.

˜ h be a function in [ H 1 (Ω)]2 . For each T ∈ Th define Theorem 5.6. Let ϕ tˆ T := σ h, T − (2ν th, T − p h, T I), and let σˆ T ∈ H(div, T ) be the unique solution of the local problem (5.5). Then, there exists C > 0, independent of h, such that

  (t, φ, p , σ , u, γ , ξ, α ) − (th , φh , ph , σ h , uh , γ , ξh , αh )  C h H

 T ∈Th

2 θ˜ T + R 2Σ

1/2 ,

(5.40)

40

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

where

  2  2 θ˜ T : = tˆ T 2[ L 2 (T )]2×2 + tr(th,T ) L 2 (T ) + σˆ T 2H(div,T ) + f + div σ h,T 2[ L 2 (T )]2 + σ h,T − σ hT,T [ L 2 (T )]2×2  2       1   + log 1 − C h (Γ ) h T Wφh + . I + K (σ h n)  2

(5.41)

[ L 2 (∂ T ∩Γ )]2

Proof. This follows from Theorem 5.2 and Lemma 5.5 by noting that h j is smaller than h T , where T ∈ Th is the triangle having Γ j as a boundary edge. 2

˜ h , and because of (5.32) (cf. Lemma 5.4), we take the idea of enforcing the quasiNow, concerning the choice of ϕ ˜ h |∂ T efficiency of θ T so that it becomes closer to full efficiency. Thus, the idea described above requires that the traces ϕ have to be as close as possible to the exact traces u|∂ T for all T ∈ Th . Since the exact solution u is unknown, the above ˜ h . Indeed, we first criterion should be understood in an empirical sense. According to this, we set an appropriate choice of ϕ ¯ 2 as the unique function satisfying the following conditions: define ϕ h ∈ [C (Ω)] 1. 2. 3. 4.

ϕ h |T ∈ [P1 ( T )]2 for all T ∈ Th ; ϕ h (x) = g (x) for each vertex x of Th lying on Γ0 ; ϕ h (x) = − V(σ h n)(x) + ( 12 I + K)ϕ h (x) + αh for each vertex x of Th lying on Γ ; ϕ h (x) is the weighted average of the constant values of uh on all the triangles T ∈ Th to which x belongs for any

interior vertex x of Th . The weighting is according to the relative area of each triangle. It follows that

ϕ h ∈ [ H 1 (Ω)]2 . Hence, we just put

ϕ˜ h |∂ T = ϕ h |∂ T ∀ T ∈ Th . Furthermore, we can also use ϕ h to deduce a simpler fully local a-posteriori error estimate, which does not require the explicit solution of the local problem (5.5). More precisely, we have the following result. Theorem 5.7. Let ϕ h be as indicated above. For each T ∈ Th we let {xk, T }k3=1 be its vertices and define

tˆ T := σ h, T − (2ν th, T − p h, T I). Then, there exists C > 0, independent of h, such that

 1/2   (t, φ, p , σ , γ , u, ξ, α ) − (th , φh , ph , σ h , uh , γ , ξh , αh )  C ˆ 2 + R2 θ , h T Σ H

(5.42)

T ∈Th

where

 2 2 θˆ T : = tˆ T 2[ L 2 (T )]2×2 + tr(th,T ) L 2 (T ) + th,T + γ h,T − ∇ ϕ h,T 2[ L 2 (T )]2×2 3     uh,T − ϕ (xk,T )2 2 + σ h,T − σ hT,T [ L 2 (T )]2×2 + h2T h, T R k =1

+ h2T |ξh |2

2 h, T [ L 2 ( T )]2

+ f + div σ + g − ϕ h 2[ H 1/2 (∂ T ∩Γ )]2 0 2      1  + I + K φh − αh   ϕ h + V(σ h n) − 2

[ H 1/2 (∂ T ∩Γ )]2

 2        φh + 1 I + K (σ h n) + log 1 − C h (Γ ) h T  W   2 2

[ L (∂ T ∩Γ )]2

Proof. This follows from Eq. (5.28) of Lemma 5.3 and by noting that

uh − ϕ h 2[ L 2 (T )]2  h2T since uh ∈ [P0 ( T )] and 2

3  uh − ϕ k =1

2  h, T (xk, T ) R2 ,

ϕ h,T ∈ [P1 ( T )]2 , which completes the proof. 2

.

(5.43)

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

41

Fig. 1. Square domain with hole.

Table 1 Uniform mesh, errors and convergence rates, Example 1. N tot

δt

1275 4915 19 299

0.63502 0.46658 0.26702

αt

δσ

0.228 0.408

5.91800 4.29526 2.55901

ασ

δu

0.238 0.379

0.24817 0.12021 0.06286

αu

δp

0.537 0.474

2.14648 1.50274 0.99632

αp

η

αη

0.264 0.300

4.63075 3.71925 2.80174

0.162 0.207

Table 2 Uniform mesh, error indicators, Example 1. N tot

ηt

ηtr

ησ σ

ηt γ ϕ

ηu−ϕ

η g −ϕ

ηξ

ηf

ηV

ηW

1275 4915 19 299

.4E−10 .2E−09 .3E−07

.3E−10 .8E−10 .2E−07

1.2100 2.0916 1.8452

1.2052 0.9524 0.5085

0.2674 0.1211 0.0559

0.5431 0.2762 0.1388

0.0045 0.0168 0.0291

.8E−11 .2E−10 .2E−08

0.1559 0.0782 0.0596

4.2587 2.9075 2.0395

Table 3 Uniform mesh, errors and convergence rates, Example 2. N tot

δt

1275 4915 19 299 76 483

3.14220 3.12946 3.02038 2.48746

αt

δσ

0.003 0.026 0.141

30.5092 29.9160 28.1725 22.6264

ασ

δu

0.015 0.044 0.159

0.48225 0.32952 0.24466 0.16038

αu

δp

0.282 0.218 0.307

12.2252 11.5797 10.2430 7.61432

αp

η

αη

0.040 0.090 0.215

8.08652 7.49701 7.31510 9.90645

0.056 0.018 −0.22

6. Numerical experiments First, we solve our model problem (2.1) on the geometry given in Fig. 1. The exact solution of our model problem is given by ν = 4 and

u(x) =

1

ν



∂x0

p (x) = 2∂x0

− log r +

(x0 −¯x0 )2

r2 (x0 −¯x0 )(x1 −¯x1 ) r2

x0 − x¯ 0 r2

=

2 r2

−

=

1



x − x¯ r2

ν

− 2(x0 −

x− x¯ 0 )2 4 r

x¯

 ,

4(x0 − x¯ 0 )2 r4

with r = |x − x¯ | and x¯ = (0, 0.5), such that  f ≡ 0 and g = u|Γ0. For this exact solution we  can observe that σ · n = O(|x|−2 ), |x| → ∞. Therefore we have 0 = B (0)−(Ω∪Ω ) div σ = ∂ B (0) σ n ds + Γ σ n ds → Γ σ n ds for R → ∞. Conse-



R

R

0

quently, there holds Σ = Γ σ · n ds = (0, 0), i.e. we have a decaying behavior at infinity. In Tables 1, 3 and 5 we use the following abbreviations for the different errors δt = t − th  L 2 (Ω) , δσ = σ − σ h  H (div;Ω) , δu = u − uh  L 2 (Ω) , δ p =  p − ph  L 2 (Ω) . The convergence rate is determined numerically with respect to the total number −α

of degrees of freedom N tot , e.g. δt ≈ C N tot t . In Tables 2, 4 and 6 we use the following notation for the individual components of the error indicator in Theorem 5.7

ηt2 =

  σ h − (2ν th − ph I)2 2

[ L ( T )]2×2

T

,

ηtr2 =

  tr(th )22

L (T )

T

,

42

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

Table 4 Uniform mesh, error indicators, Example 2. N tot

ηt

ηtr

ησ σ

ηt γ ϕ

ηu−ϕ

η g −ϕ

ηξ

ηf

ηV

ηW

1275 4915 19 299 76 483

.1E−09 .2E−09 .7E−07 .9E−06

.4E−10 .1E−09 .4E−07 .9E−06

2.2155 3.6794 5.0459 8.7248

4.0380 3.9600 4.0118 4.1613

1.2202 0.5802 0.2773 0.1375

5.7850 4.3005 2.7658 1.5518

0.0402 0.0880 0.0943 0.0911

.2E−10 .4E−10 .5E−08 .4E−07

0.1556 0.1781 0.1556 0.1438

3.0327 2.8486 2.0484 1.4979

Table 5 Adaptive version, ϑ = 0.9, errors and convergence rates, Example 2. N tot

δt

1971 3322 5391 8491 13 725 21 931

3.14220 3.13397 3.01961 2.51480 1.60281 0.96323

αt

δσ

0.005 0.077 0.403 0.938 1.086

30.5092 29.9975 28.3597 23.1882 15.8194 10.0172

ασ

δu

0.032 0.116 0.443 0.796 0.975

0.48225 0.37338 0.28372 0.22608 0.19275 0.18329

αu

δp

0.490 0.567 0.500 0.332 0.107

12.2252 11.6459 10.5051 8.15304 6.55123 4.52574

αp

η

αη

0.093 0.213 0.558 0.456 0.789

8.08652 7.47309 8.77922 13.4266 10.7049 6.42065

0.151 −0.33 −0.94 0.472 1.091

Table 6 Adaptive version, ϑ = 0.9, error indicators, Example 2. N tot

ηt

ηtr

ησ σ

ηt γ ϕ

ηu−ϕ

η g −ϕ

ηξ

ηf

ηV

ηW

1971 3322 5391 8491 13 725 21 931

.5E−14 .7E−14 .1E−13 .3E−12 .1E−11 .3E−11

.2E−14 .3E−14 .7E−14 .3E−13 .1E−12 .3E−12

2.2155 3.1812 6.7227 12.611 9.9840 5.3857

4.0380 3.9236 3.6884 3.1738 2.6528 2.4155

1.2202 0.6927 0.3272 0.1847 0.1391 0.1237

5.7850 4.3028 2.7707 1.5608 0.8844 0.3885

0.0402 0.0869 0.1239 0.1298 0.1195 0.1216

.3E−14 .5E−14 .2E−13 .4E−13 .1E−12 .3E−12

0.1556 0.2144 0.2591 0.2568 0.2371 0.2384

3.0327 3.3592 3.2264 2.9360 2.6466 2.4792

Fig. 2. Stokes (2d-FEM–BEM, Example 1): u − uh  L 2 (Ω) (left) and  p − p h  L 2 (Ω) (right).

ηt2γ ϕ =



th + γ h − ∇ ϕ h 2[ L 2 (T )]2×2 ,

ησ2 σ =

T

ηu2−ϕ =



uh − ϕh 2[ L 2 (T )]2 ,

η2g −ϕ =

T

e ⊂Γ0

 g − ϕh 2H 1/2 (e) ,

2      1 2    ηV = ϕ h + V(σ h n) − 2 I + K φh − αh  1/2

η2W

2      W  φh + 1 I + K (σ h n) =  2 2 e ⊂Γ

L (e )

H

,

σ h − σ hT 2[ L 2 (T )]2×2 ,

T



e ⊂Γ



ηξ2 =

, (e )

η2f =



h2T |ξh |2 ,

T



 div σ h − f2[ L 2 (T )]2 .

T

The errors for the first example, see Table 1 and Figs. 2 and 3 clearly indicate that the uniform version converges as predicated with order 0.5 with respect to N (or order 1 with respect to h ≈ N −1/2 ). The error indicator appears to converge sup-optimal with half the rate.

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

43

Fig. 3. Stokes (2d-FEM–BEM, Example 1): t − th  L 2 (Ω) (left) and σ − σ h  L 2 (Ω) (right).

Fig. 4. Adaptively refined meshes for Example 2.

In our second example we shift our quasi-singular point x¯ nearer to the boundary Γ0 , i.e. we investigate the previous example with the exact solution for x¯ = (0, 0.9). The adaptive version is computed with ϑ = 0.90 using a blue–green refinement strategy. In Table 3 we see the errors and their convergence rates for the uniform version, Table 4 contains the corresponding error indicators. The errors for the second example, see Table 3, show the influence of the quasi-singularity near to the boundary, which leads to a drastically longer pre-asymptotic behavior. The error indicator shows also clearly the pre-asymptotic behavior. The adaptive version shows also initially a long pre-asymptotic behavior with a very slow convergence, but it can clearly

44

M.A. Barrientos, M. Maischak / Applied Numerical Mathematics 63 (2013) 25–44

seen that, see Table 5, that the convergence in all norms, as well as the indicator, increases to the expected higher value of order 1. The sequence of meshes, see Fig. 4, shows refinement towards the singular point. As seen in Tables 4 and 6 initially the error indicators for the fem part dominate as expected, therefore most of the refinements are towards the boundary Γ0 . References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]

D.N. Arnold, F. Brezzi, J. Douglas, PEERS: A new mixed finite element method for plane elasticity, Jpn. J. Appl. Math. 1 (1984) 347–367. N.S. Bakhvalov, Solution of the Stokes nonstationary problems by the fictitious domain method, Russ. J. Numer. Anal. Math. Model. 10 (1995) 163–172. R.E. Bank, A. Weiser, Some a posteriori error estimators for elliptic partial differential equations, Math. Comp. 44 (1985) 283–301. W. Bao, Error bounds for the finite element approximation of the exterior Stokes equations in two dimensions, IMA J. Numer. Anal. 23 (2003) 125–148. W. Bao, H. Han, Error bounds for the finite element approximation of an incompressible material in an unbounded domain, Numer. Math. 93 (2003) 415–444. M.A. Barrientos, G.N. Gatica, N. Heuer, An a-posteriori error estimate for a linear-nonlinear transmission problem in plane elastostatics, Calcolo 39 (2002) 61–86. M.A. Barrientos, G.N. Gatica, M. Maischak, A-posteriori error estimates for linear exterior problems via mixed-FEM and DtN mappings, M2AN Math. Model. Numer. Anal. 36 (2) (2002) 241–272. M.A. Barrientos, G.N. Gatica, E. Stephan, A mixed finite element method for nonlinear elasticity: two-fold saddle point approach and a-posteriori error estimate, Numer. Math. 91 (2) (2002) 197–222. F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, Berlin, Heidelberg, New York, 1991. U. Brink, E.P. Stephan, Adaptive coupling of boundary elements and mixed finite elements for incompressible elasticity, Numer. Methods Partial Differential Equations 17 (1) (2001) 79–92. R. Bustinza, G.N. Gatica, M. González, A mixed finite element method for the generalized Stokes problem, Internat. J. Numer. Methods Fluids 49 (8) (2005) 877–903. J. Cahouet, J.P. Chabard, Some fast 3D finite element solvers for the generalized Stokes problem, Internat. J. Numer. Methods Fluids 8 (1988) 869–895. C. Calgaro, J. Laminie, On the domain decomposition method for the generalized Stokes problem with continuous pressure, Numer. Methods Partial Differential Equations 16 (1) (2000) 84–106. C. Carstensen, P. Causin, R. Sacco, A posteriori dual-mixed (hybrid) adaptive finite element error control for Lamé and Stokes equations, Numer. Math. 101 (2) (2005) 309–332. C. Carstensen, S. Funken, A-posteriori error control in low-order finite element discretisations of incompressible stationary flow, Math. Comp. 70 (236) (2001) 1353–1381. G.N. Gatica, Solvability and Galerkin approximations of a class of nonlinear operator equations, Z. Anal. Ihre Anwend. 21 (3) (2002) 761–781. G.N. Gatica, M. González, S. Meddahi, A low-order mixed finite element method for a class of quasi-Newtonian Stokes flows. I. A-priori error analysis, Comput. Methods Appl. Mech. Engrg. 193 (9–11) (2004) 881–892. G.N. Gatica, M. Gonzalez, S. Meddahi, A low-order mixed finite element method for a class of quasi-Newtonian Stokes flows. II. A posteriori error analysis, Comput. Methods Appl. Mech. Engrg. 193 (9–11) (2004) 893–911. G.N. Gatica, N. Heuer, A dual–dual formulation for the coupling of mixed-FEM and BEM in hyperelasticity, SIAM J. Numer. Anal. 38 (2) (2000) 380–400. G.N. Gatica, N. Heuer, S. Meddahi, On the numerical analysis of nonlinear two-fold saddle point problems, IMA J. Numer. Anal. 23 (2) (2003) 301–330. G.N. Gatica, N. Heuer, E.P. Stephan, An implicit-explicit residual error estimator for the coupling of dual-mixed finite elements and boundary elements in elastostatics, Math. Methods Appl. Sci. 23 (2001) 179–191. G.N. Gatica, A. Marquez, S. Meddahi, A new coupling of mixed finite element and boundary element methods for an exterior Helmholtz problem in the plane, Adv. Comput. Math. 30 (3) (2009) 281–301. G.N. Gatica, S. Meddahi, R. Oyarzúa, A conforming mixed finite-element method for the coupling of fluid flow with porous media flow, IMA J. Numer. Anal. 29 (1) (2009) 86–108. G.N. Gatica, E.P. Stephan, A mixed-FEM formulation for nonlinear incompressible elasticity in the plane, Numer. Methods Partial Differential Equations 18 (1) (2002) 105–128. V. Girault, P.A. Raviart, Finite Element Approximation of the Navier–Stokes Equations: Theory and Algorithms, Springer, Berlin, Heidelberg, New York, 1986. H. Han, W. Bao, The artificial boundary conditions for incompressible materials on an unbounded domain, Numer. Math. 77 (1997) 347–363. J.S. Howell, Dual-mixed finite element approximation of Stokes and nonlinear Stokes problems using trace-free velocity gradients, J. Comput. Appl. Math. 231 (2009) 780–792. G. Hsiao, W. Wendland, Boundary Integral Equations, Appl. Math. Sci., vol. 144, Springer-Verlag, Berlin, Heidelberg, 2008. C. Johnson, J.C. Nedelec, On the coupling of boundary integral and finite element methods, Math. Comp. 35 (1980) 1063–1079. G.M. Kolbelkov, M.A. Olshanskii, Effective preconditioning of Uzawa type schemes for a generalized Stokes problem, Numer. Math. 86 (2000) 443–470. A. Marquez, S. Meddahi, New implementation techniques for the exterior Stokes problem in the plane, J. Comput. Phys. 172 (2001) 685–703. S. Meddahi, F.-J. Sayas, A fully discrete BEM–FEM for the exterior Stokes problem in the plane, SIAM J. Numer. Anal. 37 (6) (2000) 2082–2102. B.V. Pal’tsev, On rapidly converging iterative methods with incomplete splitting of boundary conditions for a multidimensional singularly perturbed system of Stokes type, Russ. Acad. Sci. Sb. Math. 81 (1995) 487–531. A.J. Radcliffe, FEM–BEM coupling for the exterior Stokes problem with non-conforming finite elements and an application to small droplet deformation dynamics, Internat. J. Numer. Methods Fluids 66 (2011). P.A. Raviart, J.-M. Thomas, A mixed finite element method for 2nd order elliptic problems, in: B.E.A. Dold (Ed.), Lecture Notes in Math., vol. 606, Springer, Berlin, Heidelberg, New York, 1975. V. Sarin, A. Sameh, An efficient iterative method for the generalized Stokes problem, SIAM J. Sci. Comput. 19 (1) (1998) 206–226. A. Sequeira, The coupling of boundary integral and finite methods for the bidimensional exterior steady Stokes problem, Math. Methods Appl. Sci. 5 (1983) 356–375. A. Sequeira, On the computer implementation of a coupled boundary and finite element method for the bidimensional exterior steady Stokes problem, Math. Methods Appl. Sci. 8 (1986) 117–138.