Approximation of coupled Stokes–Darcy flow in an axisymmetric domain

Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

Contents lists available at SciVerse ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Approximation of coupled Stokes–Darcy flow in an axisymmetric domain V.J. Ervin Department of Mathematical Sciences, Clemson University, Clemson, SC 29634-0975, USA

a r t i c l e

i n f o

Article history: Received 12 September 2012 Received in revised form 21 January 2013 Accepted 12 February 2013 Available online 19 February 2013 AMS subject classifications: 76D07 76S05 74S05 65M15 35Q35

a b s t r a c t In this article we investigate the numerical approximation of coupled Stokes and Darcy fluid flow equations in an axisymmetric domain. The fluid flow is assumed to be axisymmetric. Rewriting the problem in cylindrical coordinates reduces the 3-D problem to a problem in 2-D. This reduction to 2-D requires the numerical analysis to be studied in suitably weighted Hilbert spaces. In this setting we show that the proposed approximation scheme has a unique solution, and derive corresponding a priori error estimate. Computations for an example with a know solution are presented which support the a priori error estimate. Computations are also given for a model of fluid flow in the eye. Ó 2013 Elsevier B.V. All rights reserved.

Keywords: Axisymmetric flow Stokes equation Darcy equation Coupled fluid flow

1. Introduction For the past several years the investigation of the numerical approximation of coupled Stokes–Darcy fluid flow problems has been an active area of research. A number of different formulations have been studied. An approach introduced and analyzed by Layton et al. [28] formulates the problem in terms of the unknown velocity–pressure variables in both the Stokes and Darcy domains. Other researchers who have used this approach include [38,8,19,17,20]. In [14] Discacciati et al. formulated the problem in terms of velocity–pressure unknowns in the Stokes domain and a pressure unknown (satisfying the Poisson equation) in the Darcy domain. Other researchers who have used this approach include [31,22,11,9]. Another approach is to use the Brinkman equations for the coupled problem [23]. In this approach the Stokes modeling equations and the Darcy equations are melded together using cutoff functions to obtain a single modeling equation throughout the coupled domain. Though computationally attractive this approach lacks rigorous mathematical analysis as to the relationship between the computed approximation and the true solution.

E-mail address: [email protected] 0045-7825/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2013.02.004

In this paper we investigate the special case of a 3-D Stokes– Darcy fluid flow problem in an axisymmetric domain, having an axisymmetric solution. Our motivation for considering this problem was to model fluid flow in the eye. Using the axisymmetric, we reformulate the problem in cylindrical coordinates, reducing the 3-D problem in ðx; y; zÞ to a 2-D problem in ðr; zÞ. Accompanying this reduction in spatial dimension is that the function space setting for the problem is now cast in weighted Sobolev spaces. The numerical analysis of the finite element approximation to the axisymmetric Stokes problem was presented in [4]. (See also [5,30,16,6].) For the axisymmetric Darcy problem the numerical analysis of the finite element approximation was recently given in [15]. Herein we combine the analysis of the axisymmetric Stokes and Darcy problems with the framework of [28] to obtain and analysis an approximation scheme for the coupled, axisymmetric Stokes–Darcy fluid flow problem. This paper is organized as follows. In Section 2 the modeling equations are presented, rewritten in cylindrical coordinates, and a corresponding weak formulation for the solution derived. Following, in Section 3 we present the framework for the finite element approximation, establish existence and uniqueness of the finite element approximation, and derive an a priori error estimate for the approximation. Two examples are given in the numerical experiments section. The first example investigates the a priori

97

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

error estimate for several different choices of approximating elements. The second example numerically investigates a model for fluid flow in the eye. 2. Modeling equations   R3 , denote the flow domain of interest. Additionally, let Let X  p denote bounded convex polygonal domains for the  Xf and X Stokes flow and Darcy flow, respectively. The interface boundary  f \ @X  p . Note that  :¼ @ X between the domains is denoted by C   f   X :¼ Xf [ Xp [ C. The outward pointing unit normal vectors to X  p are denoted n  let t1 ; t2 denote  f and n  p , respectively. On C and X linearly independent unit tangent vectors. We assume that there f n C  in , a subset of @ X  , which is separated is an inflow boundary C p n C  , and an outflow boundary C  out , a subset of @ X  , which from C  is also separated from C. See Fig. 2.1 for an illustration of the domain of the problem.  f n ðC  p n ðC  f :¼ @ X  [C  in Þ, and C  p :¼ @ X  [C  out Þ. Define C  p is governed We assume that the flow in the porous domain X by the Darcy’s equation subject to incompressibility of the fluid, a  out , and a non-penetration condition specified flow rate (N) across C  p. on C For the Stokes flow:

   u f Þ  p f I ¼ f f r  2m dð

f; in X

ð2:1Þ

f;  f ¼ 0 in X ru Z  C in

ð2:2Þ

 f ds ¼ N; f  n u

ð2:3Þ

f :  f ¼ 0 on C u 2

ð2:4Þ 3

ufx  f ¼ 4 ufy 5 ¼ ufx ex þ ufy ey þ ufz ez , for ex ; ey ; ez denoting unit where u ufz vectors in the x; y and z directions, respectively, and  u  Þ :¼ 1=2ðru  þ ðru  ÞT Þ represents the deformation tensor. In dð   the pressure, f an (2.1)–(2.4) u denotes the fluid’s velocity, p f

f

f

external forcing function, m the fluid kinematic viscosity, and N a specified inflow rate for the fluid.  p: For the porous domain X

 p; meff K 1 u p þ rpp ¼ f p in X

ð2:5Þ

 p;  p ¼ 0 in X ru Z  p ds ¼ N; p  n u

ð2:6Þ

 out C

 p:  p ¼ 0 on C p  n u

ð2:7Þ

p; p p ; f p , denote the fluid velocity, pressure and exterIn (2.5)–(2.8) u nal forcing functions, respectively. Additionally, in (2.5) meff repre the sents an effective kinematic fluid viscosity, and K permeability (symmetric, positive definite) tensor of the domain. For simplicity, we assume that there is a scalar function j such that  1 . jI ¼ meff K  the flows are coupled via the conservaAcross the interface C tion of mass and balance of normal forces. In addition, the Beavers–Joseph–Saffman (BJS) condition [3,24,39] is used for the . tangential forces boundary condition on C

   u f Þ  n f þ u p  n  p ¼ 0; pf  2m dð f  n f  n  f ¼ pp u    u  ; l ¼ 1; 2; f Þ  n  f  tl on C  f  tl ¼ al 2m dð u

ð2:10Þ

where a1 ; a2 denote friction constants. The boundary conditions (2.3) and (2.7) are commonly referred to as defective boundary conditions, as they do not uniquely define a solution to (2.1)–(2.8). In Section 2.1 we present a weak formulation of (2.1)–(2.8) and discuss the existence and uniqueness of the weak formulation. At the end of Section 2.1 we comment that, in addition to (2.3) and (2.7), the weak formulation implicitly im in and for u  f on C  p on poses additional boundary conditions for u  out . C 2.1. Function spaces and weak formulation In this section we introduce the function spaces needed to define the weak formulation for the coupled fluid flow problem described above.  :¼ H  ½0; 2pÞ  R3 be a bounded domain formed by Let H revolving the polygon H around the z-axis. For the axisymmetric formulation we introduce the following weighted function spaces and associated norms. For any real a and 1 6 p < 1, the space p a L ðHÞ is defined as the set of measurable functions w such that

kwka Lp ðHÞ ¼

Z

jwjp r a dx

1=p < 1;

H

where r ¼ rðxÞ is the radial coordinate of x, i.e. the distance of a point x in H from the symmetry axis. The subspace 1 L20 ðHÞ of 2 1 L ðHÞ denotes the functions q with weighted integral equal to zero, R H qrdx ¼ 0. We define the weighted Sobolev space 1 W l;p ðHÞ as the space of functions in 1 Lp ðHÞ such that their partial derivatives of order less that or equal to l belong to 1 Lp ðHÞ. Associated with 1 W l;p ðHÞ is the semi-norm j  j1 W l;p ðHÞ and norm k  k1 W l;p ðHÞ defined by

jwj1 W l;p ðHÞ ¼

l X p k@ kr @ lk z wk1 Lp ðHÞ k¼0

ð2:8Þ kwk1 W l;p ðHÞ ¼

l X k¼0

Fig. 2.1. Illustration of axisymmetric flow domain.

 ; ð2:9Þ on C

!1=p

;

!1=p jwjpW k;p ðHÞ 1

:

When p ¼ 2, we denote 1 W l;2 ðHÞ as 1 Hl ðHÞ. For a domain !  R3 we use the usual definitions and notation for Sobolev spaces. For vector functions v defined on ! expressed in cylindrical coordinates we use the notation v ðr; h; zÞ ¼ v r er þ v h eh þ v z ez ¼ ðv r ; v h ; v z Þ, where er ; eh ; ez denote unit vectors in the r; h, and z directions, respectively. Additionally, for v ðr; zÞ defined on H we use v ðr; zÞ ¼ ðv r ; v z Þ ¼ v r er þ v z ez . Let R/ denote a rotation with respect to / about the z-axis. A  is axisymmetric if v ¼v   R/ for all / 2 ð0; 2pÞ. A vector function v  is rotationally invariant if v  ¼ R/  v   R/ for all function v  Þ  H s ðH  Þ; s ¼ 0; 1; 2 denote those functions  s ðH / 2 ð0; 2pÞ. Let H  Þ which are axisymmetric, and H  Þ3 the space of rotation s ðH in Hs ðH  Þ3 vector functions. From [1,6,30] we have the ally invariant Hs ðH following two lemmas.

98

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

 Þ; s ¼ 0; 1, the mapping v  s ðH  ðr; h; zÞ 2 H  ðr; h; zÞ ! Lemma 1. For v s v ðr; zÞ 2 1 Hpffiffiffiffiffiffi ðHffi Þ is well defined for smooth functions and (up to a factor of 2p) is an isometry. Hence the lifting v ðr; zÞ 2 1 Hps ðffiffiffiffiffiffi HffiÞ;  Þ is also an isometry (up to a factor of 2p). v ðr; zÞ ! v ðr; h; zÞ 2 H s ðH 2

2

HÞ ¼ fv 2 1 H ðHÞ : @ r v =r 2 1 L ðHÞg:

 the mapping v  2 ðH  ðr; h; zÞ 2 H  ðr; h; zÞ ! v ðr; zÞ 2 Lemma 2. For v pffiffiffiffiffiffiÞ, ffi 2 H ð H Þ is (up to a factor of 2 p) an isometry. 1 þ  1=2 ðC  Þ and the isomThe trace space 1 H ðCÞ is defined using H  Þ and 1 H1 ðHÞ, and 1 H1=2 ðCÞ the dual space of  1 ðH etry between H 1=2 ðCÞ. 1H Let ra :¼ ½@=@r; @=@zT , and for v ¼ ðv r ; v z Þ, da ðv Þ :¼ 1=2ðra ðv Þþ ðra ðv ÞÞT Þ, and div axi ðv Þ :¼ ra  v þ v r =r. In addition, 1=2

Hðdiv axi ; HÞ :¼ For

n

v ¼ ðv r ; v z Þ 2 ð1 L2 ðHÞÞ2 : div axi v 2 1 L2 ðHÞ

v 2 Hðdiv axi ; HÞ; kvkHðdiv

axi ;HÞ

1H

1=2

ðCÞ as

hv  n; kiC :¼ hv  n; EC ki@ Xp ¼

For the case s ¼ 2 introduce 2 1 Hþ ð

  EC k is the axisymmetric restriction of E C k to @ Xp , satisfying kEC kk1 H1=2 ð@ Xp Þ 6 Ckkk1 H1=2 ðCÞ . Then, we define the operator v  n 2

o

ð2:11Þ

:

 :¼ kdiv axi ðv Þk21 L2 ðHÞ þ kv r k21 L2 ðHÞ þ

kv z k21 L2 ðHÞ Þ1=2 . Analogous to Lemma 2, we have the following rela Þ :¼ fv  Þ3 : r   v; H  2 L2 ðH tionship between Hðdiv axi ; HÞ and Hðdi  2 2   kHðdiv ;H Þ :¼ kdiv ðv  ÞkL2 ðH Þ þ kv x k2L2 ðH Þ þ v 2 L ðHÞg, where kv kv y k2L2 ðH Þ þ kv z k2L2 ðH Þ Þ1=2 .  Þ, the mapping  v; H r; v h; v  z Þ 2 Hðdi  ðr; h; zÞ ¼ ðv Lemma 3. For v   ðvffiffiffiffiffiffi r ; 0; ffi v z Þ ! ðv r ; v z Þ ¼ v ðr; zÞ 2 Hðdiv axi ; HÞ is (up to a factor of p 2p) an isometry. For the description of the fluid flow in Xf , we introduce the space 1 V 1 ðHÞ, a subset of 1 Hl ðHÞ, given by

n o 1 2 with norm 1 V ðHÞ ¼ w 2 1 H ðHÞ : w 2 1 L ðHÞ ;  1=2 kwk1 V 1 ðHÞ ¼ jwj21 H1 ðHÞ þ kwk21 L2 ðHÞ :

n o X p :¼ w 2 Hðdiv axi ; Xp Þ : w  njCp ¼ 0 ;

 Þ3 , Lemma 4. For v ðr; h; zÞ ¼ ðv r ; v h ; v z Þ 2 H 1 ðH r; v  z Þ ! ðv r ; v z Þ 2 1 V 1 ðHÞ  1 H1 ðHÞ. ðv

we

have

In order to incorporate the homogeneous boundary condition for the velocity on Cf , let

n o w 2 1 H1 ðXf Þ : w ¼ 0 on Cf ; n o ¼ w 2 1 V 1 ðXf Þ : w ¼ 0 on Cf :

Xf Þ ¼

1 1V }ð

and

1 1V }ð

Xf Þ

1 1 H} ð

The underlying space for the fluid velocity in Xp is Hðdiv axi ; Xp Þ. However a function w in this space may not have sufficient regularity for its trace to exist on @ Xp . Hence, the interpretation of the boundary condition w  njCp ¼ 0 needs to be carefully defined.  p Þ we have that v  p Þ. For  v; X  1=2 ðX  2H  2 Hðdi  n For v v 2 Hðdiv axi ; Xp Þ; k 2 1 H1=2 ðCÞ, following Galvis and Sarkis [19] (see also [17]), we define the operator v  np 2 1 H1=2 ðCÞ via an  1=2 ðCÞ the axisymmetextension operator EC k. Specifically, for  k2H  , define E    ric lifting of k from C to C k :¼ c u 0 , where c0 is the trace C  p Þ to H  p Þ, and u  p Þ is the weak  1 ðX  1=2 ð@ X  1 ðX  2H operator from H solution of

ð2:12Þ  p n ðC  [C  out Þ: on @ X

ð2:13Þ

ð2:15Þ

 1=2 ; k  kMp ¼ k  k1 L2 ðXp Þ : and kwkX p :¼ kdiv axi ðwÞk21 L2 ðXp Þ þ kwk21 L2 ðXp Þ ð2:16Þ Let

X :¼ X f  X p ;

and M :¼

  Z q 2 Mf  Mp : qrdx ¼ 0 ; X

and denote the dual space of X by X  . The axisymmetric weak formulation for (2.1)–(2.10) may be stated as: Given f 2 X  ; N 2 R, determine ðu; p; k; bÞ 2 X  M 1=2 ðCÞ  R2 such that, for all v 2 X and ðq; f; .Þ 2 M 1H 1=2 ð CÞ  R 2 , 1H

aðu; v Þ  bðv ; ðp; bÞÞ þ bI ðv ; kÞ ¼ ðf; v Þ;

1 bðu; ðq; .ÞÞ  bI ðu; fÞ ¼ .  N=ð2pÞ; 1

ð2:17Þ ð2:18Þ

where

aðu; v Þ :¼ af ðuf ; v f Þ þ ap ðup ; v p Þ;

bI ðv ; fÞ :¼

Z C

bðv ; ðq; .ÞÞ :

v f  nf fr ds þ hv p  np ; fiC ;

ð2:19Þ ð2:20Þ

and

af ðu; v Þ :¼ þ

Z Xf

Z C

 ur v r  2m da ðuÞ : da ðv Þ þ rdx r r

a1 as ðu  tÞðv  tÞrds;

ap ðu; v Þ :¼

Z Xp

ju  v rdx; Z

Z  vr q ra  v þ v  nf rds; rdx þ b r Xf Cin Z Z  vr q ra  v þ v  np r ds: rdx þ b bp ðv ; ðq; bÞÞ :¼ r Xp Cout

bf ðv ; ðq; bÞÞ :¼

For convenience of notation, let X f :¼ Xf Þ  Xf Þ, kv kX f ¼  1=2 2 2 2 , and M f :¼ 1 L ðXf Þ with k  kMf ¼ k  k1 L2 ðXf Þ . kv r k1 V 1 ðXf Þ þ jv z j1 H1 ðXf Þ

 p;  ¼ 0; in X  r  ru (  k; onC  ¼  =@np ¼ 0; ; @u u  out 0; onC

M p :¼ 1 L2 ðXp Þ:

¼ bf ðv f ; ðqf ; .1 ÞÞ þ bp ðv p ; ðqp ; .2 ÞÞ;

The relevance of the space 1 V 1 ðHÞ is apparent from the following lemma [1,6,30].

ð2:14Þ

where h; i@ Xp denotes the 1 L2 ð@ Xp Þ inner product, extended to a duality pairing. For the description of the fluid flow in Xp , let

1

1 1 H} ð

1   ki  ;  p; E  n hv C @ Xp 2p

ð2:21Þ ð2:22Þ ð2:23Þ ð2:24Þ

In (2.21) aas is the friction constant from the BJS condition. Conditions for the existence and uniqueness of a solution to (2.17), (2.18) are analogous to that for the discrete formulation (3.14), (3.15). Namely, (i) continuity of the operators að; Þ; bð; Þ, and bI ð; Þ, (ii) the coercivity of að; Þ (on an appropriate subspace), and (iii) that bð; Þ, and bI ð; Þ satisfy suitable inf-sup conditions (over appropriate subspaces). The continuity of the operators and the coercivity of að; Þ are straightforward to show. Below we establish the inf-sup conditions for bð; Þ and bI ð; Þ for the discrete approximation problem. Establishing the inf-sup conditions for the continuous setting can be done in a similar manner (see also [28,19]). Theorem 1. Given f 2 X  ; N 2 R, there exists a unique solution ðu; p; k; bÞ 2 X  M  1 H1=2 ðCÞ  R2 satisfying (2.17), (2.18).

99

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

2.2. Equivalence of the differential equations and weak formulation The weak formulation (2.17)–(2.24) is obtained by multiplying the differential equations by suitably smooth functions, integrating over the domain, and using Green’s theorem. Additionally, inte in and C  out (arising from using Green’s theorem) are regrals over C R R  f dS and b2 C v  p dS, respectively. For a   n placed by b1 C in v  n out smooth solution the steps used in deriving the weak formulation can be reversed to show that Eqs. (2.1)–(2.4), and (2.5)–(2.8) are satisfied. In addition, a smooth solution to (2.17), (2.18) satisfies the following boundary conditions (see [18,17]).  in be given by Let st on C

Remark. In the axisymmetric setting, for the construction of the R R R–T interpolant weighted integrals, i.e. @K . . . rds; K . . . rdx, are used (see [15]). Below we assume that m P 2; k P 1, and l 6 k. Note that the interfacial pressure approximating space Lh is contained in the space of the trace of the normal component of the velocities of X p;h , i.e. Lh  fv  np jC : v 2 X p;h g. Used in the analysis below are the following two interpolation properties. 1. From [4], there exists  a generalized  Clément interpolation operator IC : 1 H1 ðXi Þ ! PK2T i;h P s ðKÞ \ CðXi Þ such that

 u  Þn  f ¼ sn n  f þ st ; 2m dð  u  Þn f Þ  n  f . Then, smooth solutions to (2.5)–(2.8) where sn :¼ ð2m dð satisfy

 in : p f þ sn ¼ b1 On C  p ¼ b : On Cout : p

and st ¼ 0:

ð2:25Þ ð2:26Þ

2

kv  Ic v k1 L2 ðCÞ 6 Ch

1=2

kv k1 H1 ðXi Þ ;

i 2 ff ; pg;

s P 1:

ð3:6Þ

2. From the note above, and the definition of the R–T interpolant, IRT v ,

hv  np ; kh iC ¼ hIRT v  np ; kh iC ;

for all kh 2 Lh :

ð3:7Þ

Also used in the analysis are the discrete function space:

V h :¼ fv 2 X h : bI ðv h ; fÞ ¼ 0; for all f 2 Lh g;

3. Finite element approximation

ð3:8Þ

Z h :¼ fv 2 V h : bðv ; ðq; .ÞÞ ¼ 0; for all ðq; .Þ 2 Mh  R g: 2

In this section we discuss the finite element approximation to the coupled axisymmetric Stokes–Darcy system (2.17), (2.18). We focus our attention on conforming approximating spaces X f ;h  X f ; M f ;h  Mf ; X p;h  X p ; M p;h  M p ; Lh  1 H1=2 ðCÞ, where X f ;h ; M f ;h denote velocity and pressure spaces typically used for fluid flow approximations, and X p;h ; M p;h denote velocity and pressure spaces typically used for (mixed formulation) Darcy flow approximations. We begin by describing the finite element approximation framework used in the analysis. Let Xj  R2 ; j ¼ f ; p, be a polygonal domain and let T j;h be a triangulation of Xj . Thus, the computational domain is defined by

K 2 T f ;h [ T p;h :

X ¼ [K;

We assume that the triangulation is shape-regular and quasi-uniform, i.e. that there exist constants c1 ; c2 such that

c1 h 6 hK 6 c2 qK where hK is the diameter of triangle K, qK is the diameter of the greatest ball included in K, and h ¼ maxK2T f ;h [T p;h hK . We also assume that the triangulation on Xp induces the partition on C, which we denote T C;h . Let P k ðKÞ denote the space of polynomials on K of degree no greater than k, and RT k ðKÞ :¼ ðPk ðKÞÞ2 þ xP k ðKÞ denote the kth order Raviart–Thomas (R–T) elements [37,7]. Then we define the finite element spaces as follows.

X f ;h :¼ M f ;h

n

v 2 X f \ CðXf Þ2 : v jK 2 Pm ðKÞ; 8K 2 T f ;h

o ;

n o :¼ q 2 M f \ CðXf Þ : qjK 2 Pm1 ðKÞ; 8K 2 T f ;h ;

ð3:2Þ ð3:3Þ

 M p;h :¼ q 2 Mf : qjK 2 P k ðKÞ; 8K 2 T p;h ;

ð3:4Þ

 Lh :¼ f 2 CðCÞ : fjK 2 Pl ðKÞ; 8K 2 T C;h :

ð3:5Þ

v

Let

X 0f ;h :¼

n

o

v 2 X f ;h : vj@X nC ¼ 0 f

:

We have the following lemma. Lemma 5. There exists C f ;h > 0, such that

R

inf

0–qh 2Mf ;h

sup v

0 h 2X f ;h

Xf

qh div axi ðv h Þrdx kqh kMf kv h kX f

P C f ;h :

ð3:10Þ

Proof. For the case of the pressure space having mean value equal to zero the inf-sup condition (3.10) is established in [29]. As commented in [28], one can extend the inf-sup condition to the above pressure space via a local projector operator argument (see [7], Section VI.4). h For ðX p;h ; M p;h Þ Raviart–Thomas approximation spaces for the velocity and pressure, unlike in the Cartesian setting, ap ð; Þ is not coercive, with respect to the Hðdiv ; Xp Þ norm, on

Z p;h :¼ fv 2 X p;h :

Z

qdiv axi ðv Þr dx ¼ 0;

8q 2 Mp;h g:

ð3:11Þ

Xp

To compensate for this we add the term

c

Z

div axi ðuÞdiv axi ðv Þr dx

ð3:12Þ

Xp

to ap ðu; v Þ, where c > 0 is a fixed constant, and define

ap;c ðu; v Þ :¼ ap ðu; v Þ þ c

Z

div axi ðuÞdiv axi ðv Þr dx:

ð3:13Þ

Xp

2 RT k ðKÞ; 8K 2 T p;h ;

X p;h :¼



ð3:1Þ

ð3:9Þ

The spaces ðX f ;h ; Mf ;h Þ represent the Taylor–Hood pair of approximation spaces. The analysis below also holds for ðX p;h ; Mp;h Þ corresponding to the Brezzi–Douglas–Marini (BDM) approximating finite element spaces. Analogous to the continuous formulation , R we let X h :¼ X f ;h  X p;h ; and Mh :¼ q 2 Mf ;h  M p;h : X qrdx ¼ 0 .

In the approximation of Stokes and Navier–Stokes fluid flow problems in the Cartesian setting, the addition of the analogous term to (3.12) has received considerable attention recently as a means of improving the pointwise mass conservation of the approximation (see [34,35,33,27,10]). In case an axisymmetric source term, g, is modeled in the por p , i.e. in place of (2.6) we have the equation ous domain X  p , together with adding (3.12) to ap ðu; v Þ the term p ¼ g  in X ru

c

R

Xp

gdiv axi ðv Þr dx would be added to the right hand side of (3.14).

Discrete approximation problem: Given f 2 X  ; N 2 R, determine ðuh ; ph ; kh ; bÞ 2 ðX h  M h  Lh  R2 Þsuch that, for all v 2 X h and ðq; f; .Þ 2 M h  Lh  R2 ,

100

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

ac ðuh ; v Þ  bðv ; ðph ; bÞÞ þ bI ðv ; kh Þ ¼ ðf; v Þ; bðuh ; ðq; .ÞÞ  bI ðuh ; fÞ ¼ . 



1

ð3:14Þ

 Þ3  H1 ðX  Þ3 follows from the  1 ðX f 2 H Remark. The fact that v f f 1  3 uniqueness of the solution in H ðXf Þ .

ð3:15Þ

v p 2 H 1 ðX p Þ3 is constructed in an analogous manner to v f .  f ðx; y; zÞ; v p ðr; zÞ ¼ v  p ðx; y; zÞ and v f ;h ðr; zÞ ¼ Let v f ðr; zÞ ¼ v IC ðv f ðr; zÞÞ, v p;h ðr; zÞÞ ¼ IRT ðv p ðr; zÞ.



N=ð2pÞ;

1

where ac ðu; v Þ :¼ af ðuf ; v f Þ þ ap;c ðup ; v p Þ. A necessary condition for the existence and uniqueness of solutions to (3.14), (3.15) is that, for the discrete spaces X h ; M h , and Lh , the following two inf-sup conditions are satisfied. There exists constants C bh ; C X Ch > 0 .

inf

sup

ð0;0Þ–ðqh ;bÞ2M h R2 v h 2V h

inf sup

0–kh 2Lh u 2X h h

bðv h ; ðqh ; bÞÞ P C bh ; kv h kX kðq; bÞkMR2

bI ðuh ; kh Þ P C X Ch : kuh kX kkh k1 H1=2 ðCÞ

inf sup

Cin

wf ;h  nf rds þ b2

R

Cout

wp;h  np r ds

kwh kX kbkR2

0–b2R2 wh 2V h

P C RXh :

Z

¼ 0 on C;

Cin

and

vp;h  np ¼ 0

on C;

Cout

vp;h  np r ds P C 2 kv p;h kX

: ð3:20Þ

ci si ðxÞ 6 r  r0;i 6 si ðxÞ:

ð3:21Þ

8 2 > s ðxÞ; x 2 Cin ; 0 6 sin 6 jC2in j ; > < jCin j in 2 /in ðxÞ ¼ ðjCin j  sin ðxÞÞ; x 2 Cin ; jC2in j < sin ðxÞ 6 jCin j; jCin j > > : 0; otherwise:

If nf is discontinuous on Cin , e.g. Cin has corners, then /in can be constructed on any subset of Cin on which nf is continuous. Let a 2 1 H1=2 ð@ Xf Þ be given by

ð3:22Þ

ðx; y; zÞ denote the axisymmetric extension of a from @ Xf to and a  f . Introduce g 2 L2 ðX  f Þ as @X

gðx; y; zÞ ¼

 f j1=2 jX

Z  @X

inf ð0;0Þ–ðqh ;bÞ2Mh

R2

f;  on @ X ¼a   kH1=2 ð@X Þ 6 C: 6 C kgkL2 ðX Þ þ ka

v f

f

f

1=2

kv f k1 H1 ðXf Þ P C P C 1 kv f ;h kX f :

bðv h ; ðqh ; bÞÞ P C bh : k v vh 2V h h kX kðq; bÞkMR2 sup

ð3:23Þ

Proof. Let ðph ; bÞ 2 M h  R2 . We establish (3.23) via the following four steps. ^ h 2 V h to take care of b (i.e. the flow constraint.) Step 1. Construct u (Use Lemma 6). R ~ f ;h 2 X 0h such that bf ðu ~ f ;h ; ðpf ;h ; 0ÞÞ ¼ X pf ;h Step 2. Construct u f ^ f ;h ÞÞr dx. ðpf ;h  div axi ðu ~ ext 2 1 H1 ðXp Þ such that u ~ ext jC ¼ u ~ f ;h jC , and Step 3. Construct u I I R ~ ext ; ðpp;h ; 0ÞÞ ¼ X pp;h ðpp;h  div axi ðu ^ p;h ÞÞr dx . bp ðu p ~ ext , uh ¼ ðu ~ p;h :¼ IRT u ~ f ;h ; u ~ p;h Þ þ u ^ h , and verify (3.23) is Step 4. Let u satisfied.

^ h kX 6 CkbkR2 and ku Z Z ^ f ;h  nf r ds þ b1 u Cin

Cout

^ p;h  np r ds P C RXh kbk2R2 : ð3:24Þ b2 u

~ f ;h 2 X 0f ;h ; p ~f ;h 2 M f ;h satisfying Step 2. Consider u

8v 2 X 0f ;h ;

~ f ;h ; v Þ  bf ðv ; ðp ~ f ðu ~f ;h ; 0ÞÞ ¼ 0; a

~ f ;h ; ðq; 0ÞÞ ¼ ðq; pf ;h  div a ðu ^ f ;h ÞÞ; bf ðu

R

ur v r r r

8q 2 Mf ;h ;

ð3:25Þ ð3:26Þ



~ f ;h k2X 6 a ~ f ;h ; u ~ f ;h Þ ¼ bf ðu ~ f ;h ; ðp ~f ðu ~f ;h ; 0ÞÞ cku f Z ^ f ;h ÞÞrdx ~f ;h ðpf ;h  div axi ðu ¼ p Xf

  ^ f ;h Þk L2 ðX Þ ~f ;h kM kpf ;h kM þ kdiv axi ðu 6 kp f f 1 f   ^ ~ 6 kpf ;h kMf kpf ;h kMf þ Ckuf ;h kX f   ~f ;h kM kpf ;h kM þ kbkR2 : 6 Ckp f f

ð3:27Þ

From the inf-sup condition (3.10),

~f ;h kM 6 sup C f ;h kp f

f

 f kH1 ðX Þ and kv f

P C  Ch

R

n  f ds: a

 f Þk/ k 1=2 kH1=2 ð@ X Þ 6 pdiamðX  Note that ka  Þ 6 in H ð@ Xf Þ 6 C, and kg kL2 ðX f f  f Þj R a  nf dsj 6 C. pdiamðX @ Xf 3 1    f 2 H ðXf Þ such that From [19] we have that there exists v

f; r  v f ¼ g in X

ðv f  v f ;h Þ  nf r ds

Cin

~f ðu; v Þ :¼ X ra u : ra v þ where a r dx . The existence and f ~ f ;h and p ~f ;h 2 Mf ;h follows from the inf-sup condition uniqueness of u ~f ð; Þ. Next, note that (3.10) and the coercivity of a

Define /in : @ Xf ! R by

1

Cin

^ h 2 V h such that Step 1. From Lemma 6, there exists u p

 in and C  out that lies We assume that there is at most one point on C on the symmetry axis (r ¼ 0), and that such a point (if it exists) is an endpoint. For i 2 fin; outg, let si ðxÞ denote an arclength parameter on Ci . We have that there exists constants r 0;i ; ci , such that for x 2 Ci ,

aðxÞ ¼ /in ðxÞnf ;

Z

ð3:19Þ

f

Z

v f  nf r ds 

Lemma 7. For h sufficiently small, there exists C bh > 0 such that

ð3:18Þ

v f ;h  nf r ds P C 1 kv f ;h kX ;

Cin

Z

ð3:17Þ

Proof. Note that the inf-sup condition (3.18) is equivalent to showing that there exists non-zero v f ;h 2 X f ;h and v p;h 2 X p;h , and constants C 1 ; C 2 > 0, such that

v f ;h  nf

vf ;h  nf r ds ¼

Inequality (3.20) is analogously established. h

Lemma 6. There exists C RXh > 0 such that for h sufficiently small

R

Z

ð3:16Þ

Note that the inf-sup conditions (3.16), (3.17) differ from those in the Cartesian formulations in the operators involved, and the functions spaces and norms used. The following lemma is helpful in establishing (3.16).

b1

Now;

v 2X 0f ;h

¼ sup

Xf

~f ;h div axi ðv Þr dx p kv kX f

R  Xf

~ f ;h : ra v þ ra u

v 2X 0f ;h



r dx

kv kXf

v2X 0f ;h

6 sup

~ f ;h;r v u r r r

~ f ;h kX kv kX ku f f kv kX f

~ f ;h kX : ¼ ku f

ð3:28Þ

101

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

~ p;h þ u ^ p;h Þ; ðpp;h ; b2 ÞÞ bp ððu Z Z ~ p;h þ u ^ p;h Þr dx þ b2 ¼ pp;h div axi ðu

Combining (3.27) and (3.28) we obtain

~ f ;h kX ku f

  6 C kpf ;h kMf þ kbkR2 :

Step 3. Let gðxÞ ¼



ð3:29Þ

 p. of g from @ Xp to @ X  p Þ3 . Also,  1=2 ð@ X ~ f ;h ¼ 0, then g 2H As, for x 2 C; limx!C\@ Xp u

~ k 1=2  6 Cku ~ k 1  6 Cku kH 1=2 ð@ X Þ3 6 Cku ~ f ;h kX kg f ;h H ðCÞ f ;h H ðXf Þ p f

Cout

ð3:30Þ

 p Þ3 denote an extension of g  1 ðX  such that Let  z2H

j@ X ; zj@ X ¼ g p p

kH1=2 ð@ X Þ : and kzkH1 ðX p Þ 6 Ckg p

ð3:31Þ

 p Þ3 ; t 2 L2 ðX  p Þ satisfy  1 ðX  2H Let w 0 0

Z

p X

Z

p X

p   : rv  dX rw p ¼ r  wd  X q

Z

p X

Z

 p ¼ 0; 8v  p Þ3 tr  v  2 H10 ðX  dX

Z ^ f ;h  nf rds bðuh ;ðph ;bÞÞ ¼ kph k2M þ b1 u Cin Z ^ p;h  np rds P kph k2M þ C RXh kbk2R2 : ð3:38Þ þ b2 u

As kuh kX 6 Cðkph kM þ kbkR2 Þ, the stated result follows. h

ð3:32Þ ð3:33Þ

Lemma 8. There exists C X Ch > 0 such that for h sufficiently small

inf sup

 p Þ3    p ÞÞ  ðH1 ðX  p Þ3   1 ðX  tÞ 2 ðH Remark. The fact that ðw; L20 ðX 0 0 2  L ðXp Þ, follows from the uniqueness of the solution of (3.32), and 0

3

 p Þ  L2 ðX  p Þ. (3.33) in ðH10 ðX 0 From the Brezzi theory for mixed variational problems [7], we have the bound

  ^ k 2  þ kr  zk 2   L2 ðX Þ 6 C kp p;h kL2 ðX Þ þ kr  u krwk h L ðXp Þ L ðXp Þ p p   6 C kpp;h kMp þ kbkR2 þ kpf ;h kMf :

ð3:34Þ

~ ext ; u ~ p;h :¼ IRT u ~ h ¼ ðu ~ f ;h ; u ~ p;h Þ and uh ¼ u ~h þ u ^ h . Next Step 4. Let u ~ h 2 V h , i.e. bI ðu ~ h ; kh Þ ¼ 0, and hence that uh 2 V h . we show that u

Z

Proof. The analogous continuous inf-sup condition is established by, given k 2 1 H1=2 ðCÞ, constructing a suitable v ¼ ð0; v p Þ 2 X f  X p such that bI ðv ; kÞ P Ckv kX kkk1 H1=2 ðCÞ . As, for kh 2 Lh , and the Raviart–Thomas interpolant IRT v p , we have hv p  np ; kh iC ¼ hIRT v p  np ; kh iC , then for v h ¼ ð0; IRT v p Þ, (3.39) follows. h

Corollary 1. There exists a constant C c > 0 such that inf

sup

ð0;0;0Þ–ðqh ;fh ;.h Þ2M h Lh R2 v h 2X h

bðv h ; ðqh ; .h ÞÞ  bI ðv h ; fh Þ P Cc : ðkqh kM þ kfh k1 H1=2 ðCÞ þ k.h kR2 Þkv h kX

ð3:40Þ

Proof. The combined inf-sup condition (3.40) follows from the individual inf-sup conditions (3.23) and (3.39) (see [21,17]). h

Theorem 2. Given f 2 X  ; N 2 R, for c > 0 there exists a unique solution ðuh ; ph ; kh ; bÞ 2 ðX h  M h  Lh  R2 Þ satisfying (3.14), (3.15).

C

C

ð3:39Þ

With the inf-sup conditions given in (3.16) and (3.17) now established, we have the following.

~ p;h  np ; EC kh i@ Xp ~ f ;h  nf kh rds þ hu u ZC Z ~ f ;h  nf kh rds þ IRT u ~ ext  np kh r ds ¼ u ~ ext  np ¼ 0 on @ Xp n CÞ ðas IRT u Z Z ~ f ;h  nf kh r ds þ ðw þ zÞ  np kh r ds ¼ u ZC ZC ~ f ;h  nf kh rds þ u ~ f ;h  np kh r ds ¼ 0: u ¼

bI ðuh ; kh Þ P C X Ch : kuh kX kkh k1 H1=2 ðCÞ

From (3.23) and (3.39) we obtain the following ‘‘combined’’ infsup condition.

 p to Xp . From ~ ext denote the reduction of ðw  þ Let u zÞ from X ~ ext 2 1 V 1 ðXp Þ  1 H1 ðXp Þ with ku ~ ext kX 6 (3.30), (3.31) and (3.34), u p ~ ext k V 1 ðX Þ H1 ðX Þ 6 Cðkpf ;h kM þ kpp;h kMp þ kbkR2 Þ. Cku p p f 1 1

C

Cout

ð3:37Þ Thus

0–kh 2Lh u 2X h h

~ h ; kh Þ ¼ bI ðu

Xp

Cout

^  r  zÞdX  p ; 8q  p Þ:  ðp  2 L20 ðX p;h  r  u q h

p X

~ p;h þ u ^ p;h Þ  np r ds ðu

Cout

Xp

x 2 @ Xp n C  denote the lifting , and g x2C

6 Cðkpf ;h kMf þ kbkR2 Þ:

Z

^ p;h Þr dx pp;h div axi ðIRT ðw þ zÞÞr dx þ pp;h div axi ðu Xp Z ^ p;h  np r ds u þ b2 Cout Z Z ^ p;h Þr dx ¼ pp;h div axi ðw þ zÞr dx þ pp;h div axi ðu Xp Z Xp Z Z ^ p;h  np r ds ¼ ^ p;h  np r ds: p2p;h r dx þ b2 u u þ b2 ¼

0; ~ f ;h ; u

ZXp

ð3:35Þ

Proof. The proof follows from the continuity of the operators ac ð; Þ; bð; Þ, and bI ð; Þ, the coercivity of ac ð; Þ on Z h  Z h , and the inf-sup conditions (3.16) and (3.17). h

C

Now,

3.1. A priori error estimate sup vh 2X h

bðv h ; ðph ; bÞÞ bðuh ; ðph ; bÞÞ P kv h kX kuh kX ¼

~ f ;h þ u ^ f ;h Þ; ðpf ;h ; b1 ÞÞ þ bp ððu ~ p;h þ u ^ p;h Þ; ðpp;h ; b2 ÞÞ bf ððu : kuh kX

~ f ;h þ u ^ f ;h Þ; ðpf ;h ; b1 ÞÞ bf ððu Z Z ~ f ;h þ u ^ f ;h Þr dx þ b1 ¼ pf ;h div axi ðu Xf

¼

Z

Xf

p2f ;h r dx þ b1

Z Cin

~ f ;h þ u ^ f ;h Þ  nf r ds ðu

Cin

^ f ;h  nf r ds: u

ð3:36Þ

In this section we present the a priori error estimate for the approximation ðuh ; ph Þ. Theorem 3. For ðu; p; k; bÞ satisfying (2.17), (2.18) and ðuh ; ph ; kh ; bh Þ satisfying (3.14), (3.15), and h sufficiently small, there exists a constant C > 0 such that

ku  uh kX þ kp  ph kM þ kb  bh kR2 þ kk  kh k1 H1=2 ðCÞ   6 C inf vh 2X h ku  v h kX þ inf qh 2Mh kp  qh kM þ inf fh 2Lh kk  fh k1 H1=2 ðCÞ : ð3:41Þ

102

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

Proof. Note that the continuous solution u determined by (2.17), (2.18) satisfies

c

0.5 0.55 0.45

Z

0.5

div axi ðuÞdiv axi ðv Þr dx ¼ 0;

8v 2 X:

0.4 0.45

Xp

Hence, the continuous problem (2.17), (2.18) can be equivalently stated with aðu; v Þ replaced by ac ðu; v Þ. Then (3.14), (3.15) is simply a Galerkin discretization of a mixed variational problem in the space X  M  R2  1 H1=2 ðCÞ using the spaces X h  Mh  R2  Lh . The stated error estimate then holds from the Brezzi theory for mixed variational problems [7]. h

0.35 0.4 0.3 0.35 0.25 0.3 0.2 0.25 0.15

Remark. For ðu; pÞ sufficiently smooth, X f ;h ; M f ;h ; X p;h ; M p;h given by (3.1)–(3.4), we have from [4,15,29] that m

m

inf kuf  v h kXf 6 Ch ;

inf kpf  qh k1 L2 ðXf Þ 6 Ch ;

v h 2X f ;h

inf kup  v h kX p 6 Ch

kþ1

v h 2X p;h

ð3:42Þ

qh 2Mf ;h

inf kpp  qh k1 L2 ðXp Þ 6 Ch

;

kþ1

qh 2Mp;h

ð3:43Þ

0.15

0.05 0

:

0.2

0.1

0.1 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

Fig. 4.2. Computed Stokes velocity flow field for Example 1.

Remark. For a sufficiently smooth function k, the interpolation results can be extended to obtain

0 0.55

inf kk  fh k1 Hs ðCÞ 6 Ch

lþ1s

fh 2Lh

;

s ¼ 0; 1;

ð3:44Þ

−0.05 0.5 −0.1

and then by operator interpolation to yield

inf kk  fh k1 H1=2 ðCÞ 6 Ch

lþ1=2

fh 2Lh

0.45 −0.15

ð3:45Þ

:

0.4 −0.2 0.35 −0.25

4. Numerical experiments

0.3 −0.3

In this section we numerically investigate the approximation of two, cylindrically symmetric, coupled Stokes–Darcy flow problems. The first experiment is performed on an example with a known solution. Rates of convergence of the approximation to the known solution are computed for several different choices of approximating elements and compared with those predicted by Theorem 3. The second example we investigate is that of fluid flow through the eye. For this example we compare the flow profiles obtained assuming a parabolic inflow and outflow profile with that from using the defective boundary condition discussed above. In describing the approximation spaces/approximating elements below, we use P k to denote the space of polynomials of

0.25 −0.35

0.2

−0.4

0.15

−0.45 −0.5 0

0.1 0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

Fig. 4.3. Computed Darcy velocity flow field for Example 1.

0.5 0.4 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Fig. 4.1. Computational mesh corresponding to h ¼ 1=4.

0.45

0.5

0.5

Fig. 4.4. Illustration of the anatomy of an eye [26].

103

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

degree 6 k on each triangle, which are continuous over the domain. The notation discPk refers to the approximation spaces/ approximating elements which are polynomials of degree 6 k on each triangle, and are not required to be continuous over the domain. For all the computations presented below the value used for c in (3.13) was c ¼ 1:0. 4.1. Example 1 For this example we take X ¼ ð0; 1=2Þ  ð1=2; 1=2Þ; Xf ¼ ð0; 1=2Þ  ð0; 1=2Þ, Xp ¼ ð0;1=2Þ  ð1=2; 0Þ, and C ¼ ð0;1=2Þ  f0g. For the fluid velocity in Xf and Xp we use

"

uf ðr; zÞ ¼ up ðr; zÞ ¼

#

r cosðprÞ sinðpzÞ  p2 cosðprÞ cosðpzÞ þ r sinðprÞ cosðpzÞ

; ð4:1Þ

and for the pressure in Xf and Xp

pf ¼ pp ¼ sinðpzÞð cosðprÞ þ 2pr sinðprÞÞ þ 4re4r cosðpzÞ 

2

p

ð1  5e2 Þ:

ð4:2Þ

In addition we use m ¼ meff ¼ 1; j ¼ 1, and aas ¼ 1. Computations were performed on a series of meshes. Illustrated in Fig. 4.1 is the computational mesh corresponding to mesh parameter h ¼ 1=4. The computed flow field and contour lines for the magnitude of the velocity on the mesh h ¼ 1=8, using Taylor– Hood P 2  P1 elements for approximating ðuf ; pf Þ, Raviart–Thomas RT 1  discP1 elements for approximating ðup ; pp Þ, and P1 elements for approximating the interfacial pressure k, are given in Figs. 4.2 and 4.3. Presented in Table 4.1 are computations obtained using the mini-element approximation pair, ðP1 þ BubbleÞ  P1 for ðuf ;h ; pf ;h Þ, Raviart–Thomas RT 1  discP1 for ðup;h ; pp;h Þ, and P 1

Table 4.1 Example 1 using ðP1 þ BubbleÞ  P 1 ; RT 1  discP 1 , and P1 approximation elements. Stokes flow approximation errors h

kuf  uf ;h kX f

Cvg. rate

kuf  uf ;h k1 L2 ðXf Þ

Cvg. rate

kpf  pf ;h k1 L2 ðXf Þ

Cvg. rate

1/2 1/4 1/8 1/12 1/16

1.541E01 7.513E02 3.708E02 2.460E02 1.841E02

1.04 1.02 1.01 1.01

4.705E03 1.230E03 3.106E04 1.384E04 7.791E05

1.94 1.98 1.99 2.00

9.894E02 3.141E02 1.059E02 5.619E03 3.593E03

1.66 1.57 1.56 1.55

Darcy flow approximation errors h kup  up;h kX p

Cvg. rate

kup  up;h k1 L2 ðXp Þ

Cvg. rate

kpp  pp;h k1 L2 ðXp Þ

Cvg. rate

1/2 1/4 1/8 1/12 1/16

1.92 1.92 2.06 2.12

7.424E03 2.030E03 5.405E04 2.314E04 1.232E04

1.87 1.91 2.09 2.19

2.989E02 6.569E03 1.594E03 7.023E04 3.932E04

2.19 2.04 2.02 2.02

Cvg. rate

kk  kh k1 L2 ðCÞ

Cvg. rate

0.91 0.93 0.97 0.98

4.107E02 8.929E03 2.161E03 9.513E04 5.325E04

2.20 2.05 2.02 2.02

8.967E03 2.365E03 6.230E04 2.703E04 1.468E04

Interfacial pressure approximation errors h kk  kh k1 H1 ðCÞ 1/2 1/4 1/8 1/12 1/16

1.508E01 8.033E02 4.213E02 2.845E02 2.145E02

Table 4.2 Example 1 using P 2  P1 ; RT 1  discP 1 , and P1 approximation elements. Stokes flow approximation errors h

kuf  uf ;h kX f

Cvg. rate

kuf  uf ;h k1 L2 ðXf Þ

Cvg. rate

kpf  pf ;h k1 L2 ðXf Þ

Cvg. rate

1/2 1/4 1/8 1/12 1/16

1.761E02 4.452E03 1.115E03 4.956E04 2.788E04

1.98 2.00 2.00 2.00

5.184E04 6.693E05 8.390E06 2.486E06 1.048E06

2.95 3.00 3.00 3.00

1.450E02 3.265E03 7.836E04 3.454E04 1.937E04

2.15 2.06 2.02 2.01

Cvg. rate

kup  up;h k1 L2 ðXp Þ

Cvg. rate

kpp  pp;h k1 L2 ðXp Þ

Cvg. rate

1.94 1.93 2.06 2.13

7.415E03 2.007E03 5.331E04 2.277E04 1.210E04

1.89 1.91 2.10 2.20

7.880E03 2.032E03 5.140E04 2.291E04 1.290E04

1.96 1.98 1.99 2.00

Interfacial pressure approximation errors h kk  kh k1 H1 ðCÞ

Cvg. rate

kk  kh k1 L2 ðCÞ

Cvg. rate

1/2 1/4 1/8 1/12 1/16

0.86 0.93 0.97 0.98

4.978E03 1.286E03 3.390E04 1.529E04 8.652E05

1.95 1.92 1.96 1.98

Darcy flow approximation errors h kup  up;h kX p 1/2 1/4 1/8 1/12 1/16

8.966E03 2.344E03 6.164E04 2.671E04 1.449E04

1.451E01 7.995E02 4.210E02 2.844E02 2.145E02

104

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

Table 4.3 Example 1 using P2  P1 ; RT 2  discP 2 , and P1 approximation elements. Stokes flow approximation errors h

kuf  uf ;h kX f

1/2 1/4 1/8 1/12 1/16

1.761E02 4.452E03 1.115E03 4.956E04 2.788E04

Cvg. rate

kuf  uf ;h k1 L2 ðXf Þ

Cvg. rate

kpf  pf ;h k1 L2 ðXf Þ

Cvg. rate

1.98 2.00 2.00 2.00

5.185E04 6.692E05 8.389E06 2.485E06 1.048E06

2.95 3.00 3.00 3.00

1.451E02 3.264E03 7.836E04 3.454E04 1.937E04

2.15 2.06 2.02 2.01

Darcy flow approximation errors h kup  up;h kX p

Cvg. rate

kup  up;h k1 L2 ðXp Þ

Cvg. rate

kpp  pp;h k1 L2 ðXp Þ

Cvg. rate

1/2 1/4 1/8 1/12 1/16

1.35 1.38 1.44 1.47

2.084E02 8.177E03 3.129E03 1.743E03 1.144E03

1.35 1.39 1.44 1.46

1.277E03 1.708E04 2.532E05 9.069E06 4.536E06

2.90 2.75 2.53 2.41

Interfacial pressure approximation errors h kk  kh k1 H1 ðCÞ

Cvg. rate

kk  kh k1 L2 ðCÞ

Cvg. rate

1/2 1/4 1/8 1/12 1/16

0.87 0.92 0.96 0.98

4.971E03 1.294E03 3.407E04 1.535E04 8.673E05

1.94 1.93 1.97 1.98

2.088E02 8.199E03 3.143E03 1.751E03 1.149E03

1.452E01 7.948E02 4.193E02 2.838E02 2.142E02

Table 4.4 Example 1 using P2-P1, RT2 – discP2, and P2 approximation elements. Stokes flow approximation errors h

kuf  uf ;h kX f

Cvg. rate

kuf  uf ;h k1 L2 ðXf Þ

Cvg. rate

kpf  pf ;h k1 L2 ðXf Þ

Cvg. rate

1/2 1/4 1/8 1/12 1/16 1/20

1.761E02 4.452E03 1.115E03 4.955E04 2.788E04 1.784E04

1.98 2.00 2.00 2.00 2.00

5.163E04 6.686E05 8.387E06 2.485E06 1.048E06 5.367E07

2.95 2.99 3.00 3.00 3.00

1.450E02 3.264E03 7.835E04 3.454E04 1.937E04 1.238E04

2.15 2.06 2.02 2.01 2.01

Darcy flow approximation errors h kup  up;h kX p

Cvg. rate

kup  up;h k1 L2 ðXp Þ

Cvg. rate

kpp  pp;h k1 L2 ðXp Þ

Cvg. rate

1/2 1/4 1/8 1/12 1/16 1/20

1.85 2.22 2.36 2.41 2.43

3.029E03 8.413E04 1.806E04 6.954E05 3.483E05 2.026E05

1.85 2.22 2.35 2.40 2.43

1.143E03 1.304E04 1.474E05 4.219E06 1.751E06 8.877E07

3.13 3.14 3.09 3.06 3.04

Interfacial pressure approximation errors h kk  kh k1 H1 ðCÞ

Cvg. rate

kk  kh k1 L2 ðCÞ

Cvg. rate

1/2 1/4 1/8 1/12 1/16 1/20

1.88 1.95 1.99 2.00 2.00

1.412E03 1.819E04 2.314E05 6.933E06 2.941E06 1.510E06

2.96 2.97 2.97 2.98 2.99

3.082E03 8.579E04 1.837E04 7.053E05 3.527E05 2.049E05

3.470E02 9.454E03 2.441E03 1.091E03 6.145E04 3.934E04

approximation elements for kh . Theorem 3 predicts (bounded by the inf vf ;h 2X f ;h kuf  v f ;h k1 H1 ðXf Þ term) Fig. 4.3

kuf  uf ;h kX f þ kpf  pf ;h k1 L2 ðXf Þ þ kup  up;h kX p 1

þ kpp  pp;h k1 L2 ðXp Þ þ kk  kh k1 H1=2 ðCÞ 6 Ch :

ð4:3Þ

RT 0  discP 0 is not an appropriate choice as an approximation pair for ðup;h ; pp;h Þ. Table 4.2 contains the computations obtained using Taylor– Hood P2  P 1 approximating elements for ðuf ;h ; pf ;h Þ, Raviart– Thomas RT 1  discP1 for ðup;h ; pp;h Þ, and P 1 approximation elements for kh . Theorem 3 predicts (bounded by the inf fh 2Lh kk  fh k1 H1=2 ðCÞ term)

The experimental convergence rates computed in Table 4.1 are consistent with (4.3). Remark. As Lh  1 H1=2 ðCÞ, for a conforming method, we require that our approximation for k; kh , be a continuous function. Assumption A1, Lh  fv  np jC : v 2 X p;h g, then implies that

kuf  uf ;h kXf þ kpf  pf ;h k1 L2 ðXf Þ þ kup  up;h kX p þ kpp  pp;h k1 L2 ðXp Þ þ kk  kh k1 H1=2 ðCÞ 6 Ch

3=2

:

ð4:4Þ

The experimental convergence rates computed in Table 4.2 are consistent with (4.4).

105

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

Table 4.5 Example 2 using P2  P 1 ; RT 1  discP 1 , and P 1 approximation elements, with a defective boundary condition.

z 50 µ m

7.2 300 µ m

Canal of Schlemm 300 µ m

(1.91 , 5.47)

100 µ m

(5.9 , 4.13) 4.5

(1.53 , 4.66)

Stokes flow norms h

kuf ;h kX f

kuf ;h k1 L2 ðXf Þ

kpf ;h k1 L2 ðXf Þ

1 1/2 1/4

1.717E01 1.725E01 1.726E01

1.616E02 1.622E02 1.623E02

3.977E+00 3.981E+00 3.981E+00

Darcy flow norms h kup;h kX p

kup;h k1 L2 ðXp Þ

kpp;h k1 L2 ðXp Þ

1 1/2 1/4

1.002E03 9.967E04 9.943E04

1.363E+01 1.353E+01 1.348E+01

(6.22 , 3.62)

r

Fig. 4.5. Model for simulating fluid flow through the eye.

1.002E03 9.967E04 9.943E04

Interfacial pressure norms h kkh k1 H1 ðCÞ

kkh k1 L2 ðCÞ

1 1/2 1/4

9.637E01 1.078E+00 1.026E+00

3.017E+01 6.345E+01 7.978E+01

7.5

Table 4.6 Example 2 using P2  P 1 ; RT 1  discP 1 , and P1 approximation elements, assuming a parabolic inflow and outflow profile.

7 6.5

Stokes flow norms

6 5.5

h

kuf ;h kX f

kuf ;h k1 L2 ðXf Þ

kpf ;h k1 L2 ðXf Þ

1 1/2 1/4

1.732E01 1.742E01 1.743E01

1.620E02 1.627E02 1.628E02

4.003E+00 4.008E+00 4.009E+00

Darcy flow norms h kup;h kX p

kup;h k1 L2 ðXp Þ

kpp;h k1 L2 ðXp Þ

1 1/2 1/4

1.111E03 1.111E03 1.111E03

1.606E+01 1.607E+01 1.607E+01

5 4.5 4 3.5

0

1

2

3

4

5

6

7

Fig. 4.6. Computed mesh for the AC corresponding to h ¼ 1.

1.111E03 1.111E03 1.111E03

Interfacial pressure norms h kkh k1 H1 ðCÞ

kkh k1 L2 ðCÞ

1 1/2 1/4

1.064E+00 1.212E+00 1.154E+00

3.395E+01 7.273E+01 9.212E+01

4.3

4.2

kuf  uf ;h kX f þ kpf  pf ;h k1 L2 ðXf Þ þ kup  up;h kX p þ kpp  pp;h k1 L2 ðXp Þ 4.1

þ kk  kh k1 H1=2 ðCÞ 6 Ch

3=2

;

ð4:5Þ

and for the results in Table 4.4, Theorem 3 predicts (bounded by the inf vf ;h 2X f ;h kuf  v f ;h k1 H1 ðXf Þ and inf qf ;h 2Mf ;h kpf  qf ;h k1 L2 ðXf Þ term)

4

3.9

kuf  uf ;h kX f þ kpf  pf ;h k1 L2 ðXf Þ þ kup  up;h kX p þ kpp  pp;h k1 L2 ðXp Þ 2

þ kk  kh k1 H1=2 ðCÞ 6 Ch :

3.8

ð4:6Þ

The experimental convergence rates computed in Tables 4.3 and 4.4 are consistent with (4.5) and (4.6), respectively.

3.7

3.6 5.9

5.95

6

6.05

6.1

6.15

6.2

6.25

6.3

6.35

6.4

Fig. 4.7. Computed mesh for the TM corresponding to h ¼ 1.

The computations in Tables 4.3 and 4.4 were obtained using Taylor–Hood P 2  P1 approximating elements for ðuf ;h ; pf ;h Þ and Raviart–Thomas RT 2  discP 2 for ðup;h ; pp;h Þ. For Table 4.3, k was approximated using P1 elements, and for Table 4.4, k was approximated using P2 elements. For the results in Table 4.3, Theorem 3 predicts (bounded by the inf fh 2Lh kk  fh k1 H1=2 ðCÞ term)

Remark. The influence of using P2 elements for kh when using RT 2  discP2 for ðup;h ; pp;h Þ is clearly demonstrated by comparing the numerical results in Tables 4.3 and 4.4. 4.2. Example 2 To demonstrate the difference between flow fields generated assuming a parabolic inflow and outflow profile with that from the defective boundary conditions given by (2.3) and (2.7), we sim-

106

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108 x 10-3

Velocity in the TM region 4.1

10

4.05

9

4

8

3.95

7

3.9

6

3.85

5

3.8

4 3

3.75

2

3.7

1 3.65 5.95

Fig. 4.8. Computed flow field in the AC, h ¼ 1=4, for a defective boundary condition.

6

6.05

6.1

6.15

6.2

6.25

6.3

Fig. 4.10. Computed flow field in the TM, h ¼ 1=4, for a defective boundary condition.

Pressure in the AC region 7.5

Pressure in the TM region 3

4.3

7 2.5 6.5

4.2

−10

2 1.5

6

1 5.5

4.1 −20 4 −30

0.5 5

0

3.9 −40

−0.5

4.5

3.8

−1

−50

4 −1.5 3.5

0

1

2

3

4

5

6

3.7 −60

7

Fig. 4.9. Computed pressure profile in the AC, h ¼ 1=4, for a defective boundary condition.

ulate fluid flow in the eye. In this simulation we assume the eye is looking straight up. Fluid in the eye is generated by the ciliary body which is located on the wall of the eye adjacent to the lens [32]. Fluid flows from the ciliary body into the Anterior Cavity (front section of the eye), AC, by passing between the lens and the iris and then flowing through the pupil, see Fig. 4.5. We assume a flow rate of 2 ll /min [25]. The fluid exits the AC through the Trabecular Meshwork, TM, located on the wall of the eye slightly above where the iris attaches to the wall. After flowing through the TM the fluids enters the Canal of Schlemm. The model geometry of the eye illustrated in Fig. 4.5 was constructed using [12,13,36]. For the model we assume: 1. The radius of the inside of the cornea is 7.2 mm. 2. The radius of the lens is 12.5 mm. 3. The distance between the lens and the cornea along the vertical axis is 2.7 mm. 4. The pupil aperture is 3 mm. 5. The lower side of the iris has the same curvature as the lens. The top side of the iris is approximated as a straight line. The

3.6 5.9

5.95

6

6.05

6.1

6.15

6.2

6.25

6.3

6.35

6.4

Fig. 4.11. Computed pressure profile in the TM, h ¼ 1=4, for a defective boundary condition.

6. 7. 8. 9. 10.

width of the iris is approximately 0.5 mm, and we assume the iris attaches to the cornea. (Physically the iris attaches to the ciliary muscles very near the cornea) The gap between the iris and the lens is 0.25 mm. The length of the interface between the AC and TM is 0.6 mm. The width of the TM at the bottom is 0.1 mm. The length of the interface of the TM with the Canal of Schlemm is 0.3 mm. A straight line connects the point at the top of the interface of the AC and TM with the point at the top of the interface of the TM with the Canal of Schlemm.

The fluid flow in the AC is modeling using the Stokes equations and that through the TM modeled using the Darcy equations. For the kinematic viscosity of the fluid we use m = 0.66 mm/s (approximately that of water), and in the TM for the effective viscosity meff ¼ 0:66 mm=s, and for the permeability j ¼ 2:0  106 mm2 [2].

107

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

x 10-3

Velocity in the TM region 4.1

4

4.05

3.5

4 3 3.95 2.5 3.9 2

3.85 3.8

1.5

3.75

1

3.7 0.5 3.65 5.95

Fig. 4.12. Computed flow field in the AC, h ¼ 1=4, for a parabolic inflow and outflow profile.

6

6.05

6.1

6.15

6.2

6.25

6.3

Fig. 4.14. Computed flow field in the TM, h ¼ 1=4, for a parabolic inflow and outflow profile.

Pressure in the TM region

Pressure in the AC region

4.3

7.5 3

−10 4.2

7 2.5 6.5

2

−20 4.1 −30

1.5

6

1 5.5

0.5

4

−40 −50

3.9

0

5

−0.5

−60 3.8 −70

4.5 −1

3.7 4

−80

−1.5

−90 3.5

−2 0

1

2

3

5

6

7

Fig. 4.13. Computed pressure profile in the AC, h ¼ 1=4, for a parabolic inflow and outflow profile.

In our modeling equations there arises from the Beavers–Joseph–Saffman boundary condition (2.10) a frictional constant, denoted by aas in (2.21). As we do not have a reference value for aas , we use the nominal value aas ¼ 1 in our computations. Computations were performed on a series of three meshes. The initial mesh (h ¼ 1) is given in Figs. 4.6 and 4.7. The two subsequent meshes (h ¼ 1=2; h ¼ 1=4) were obtained by mid edge refinement, whereby each triangle in the mesh was divided into four smaller triangles by connecting the mid edges of each side. The approximating elements used were the Taylor–Hood P 2  P 1 for the velocity–pressure approximation in the AC, Raviart–Thomas RT 1  discP1 for the velocity–pressure in the TM, and continuous, piecewise linear elements, P1 , for the interfacial pressure term along the interface between the AC and TM regions. Presented in Tables 4.5 and 4.6 are various norms of the approximation computed on the three different meshes for the case of the degenerate boundary condition. Illustrated in Figs. 4.8, 4.9, 4.10 and 4.11, 4.12, 4.13, 4.15 are the flow fields and pressure plots obtained on the finest mesh, h ¼ 1=4, for the cases of the degenerate boundary condition and an assumed

3.6 5.9

5.95

6

6.05

6.1

6.15

6.2

6.25

6.3

6.35

6.4

Fig. 4.15. Computed pressure profile in the TM, h ¼ 1=4, for a parabolic inflow and outflow profile.

parabolic inflow and outflow profile, respectively. Visually the difference in the computed solutions is most apparent in the velocity and pressure plots in the TM, i.e. by comparing Figs. 4.10, 4.11, with Figs. 4.14, 4.15. Consistent with (2.26), for simulations generated using the defective boundary condition formulation we note that the pressure along the outflow boundary is constant, and the flow much more uniform than for the case when the parabolic profile for the outflow is imposed. Fig. 4.4. Acknowledgement We gratefully acknowledge the suggestions made by the referees, which have significantly improved this paper. References [1] F. Assous, P. Ciarlet Jr., S. Labrunie, Theoretical tools to solve the axisymmetric Maxwell equations, Math. Methods Appl. Sci. 25 (1) (2002) 49–78. [2] R. Avtar, R. Srivastava, Modelling aqueous humor outflow through trabecular meshwork, Appl. Math. Comput. 189 (2007) 34–745. [3] G. Beavers, D. Joseph, Boundary conditions at a naturally impermeable wall, J. Fluid Mech. 30 (1967) 197–207.

108

V.J. Ervin / Comput. Methods Appl. Mech. Engrg. 258 (2013) 96–108

[4] Z. Belhachmi, C. Bernardi, S. Deparis, Weighted Clément operator and application to the finite element discretization of the axisymmetric Stokes problem, Numer. Math. 105 (2) (2006) 217–247. [5] Z. Belhachmi, C. Bernardi, S. Deparis, F. Hecht, A truncated Fourier/finite element discretization of the Stokes equations in an axisymmetric domain, Math. Models Methods Appl. Sci. 16 (2) (2006) 233–263. [6] C. Bernardi, M. Dauge, Y. Maday, Spectral methods for axisymmetric domains, volume 3 of Series in Applied Mathematics (Paris), Gauthier-Villars, 1999. [7] F. Brezzi, M. Fortin, Mixed and hybrid finite element methods, Springer-Verlag, New York, 1991. [8] E. Burman, P. Hansbo, A unified stabilized method for Stokes’ and Darcy’s equations, J. Comput. Appl. Math. 198 (1) (2007) 35–51. [9] Y. Cao, M. Gunzburger, X. Hu, F. Hua, X. Wang, W. Zhao, Finite element approximations for Stokes–Darcy flow with Beavers–Joseph interface conditions, SIAM J. Numer. Anal. 47 (6) (2010) 4239–4256. [10] M.A. Case, V.J. Ervin, A. Linke, L.G. Rebholz, A connection between Scott– Vogelius and grad-div stabilized Taylor–Hood FE approximations of the Navier–Stokes equations, SIAM J. Numer. Anal. 49 (2011) 1461–1481. [11] A. Çesßmeliog˘lu, B. Rivière, Primal discontinuous Galerkin methods for timedependent coupled surface and subsurface flow, J. Sci. Comput. 40 (1-3) (2009) 115–140. [12] M.W. Charles, N. Brown, Dimensions of the human eye relevant to radiation protection, Phys. Med. Biol. 20 (1975) 202–218. } ke, G.K. Kreiglstein, Morphological variability of [13] T.S. Dietlein, P.C. Jacobi, C. Lu the trabecular meshwork in glaucoma patients: implications for nonperforating glaucoma surgery, Br. J. Ophthalmol. 84 (2000) 1354–1359. [14] M. Discacciati, E. Miglio, A. Quateroni, Mathematical and numerical models for coupling surface and groundwater flows, Appl. Numer. Math. 43 (2002) 57–74. [15] V.J. Ervin, Approximation of axisymmetric Darcy flow. Technical report, Dept. Math. Sci., Clemson University, 2012. . [16] V.J. Ervin, E.W. Jenkins, Stenberg’s sufficiency condition for axisymmetric Stokes flow. Technical report, Dept. Math. Sci., Clemson University, 2011. . [17] V.J. Ervin, E.W. Jenkins, S. Sun, Coupled generalized non-linear Stokes flow with flow through a porous media, SIAM J. Numer. Anal. 47 (2009) 929–952. [18] L. Formaggia, J.-F. Gerbeau, F. Nobile, A. Quateroni, Numerical treatment of defective boundary conditions for the Navier–Stokes equations, SIAM J. Numer. Anal. 40 (2002) 376–401. [19] J. Galvis, M. Sarkis, Non-matching mortar discretization analysis for the coupling Stokes–Darcy equations, Electron. Trans. Numer. Anal. 26 (2007) 350–384. [20] G.N. Gatica, R. Oyarzúa, F.-J. Sayas, Analysis of fully-mixed finite element methods for the Stokes–Darcy coupled problem, Math. Comput. 80 (276) (2011) 1911–1948.

[21] G.N. Gatica, F.-J. Sayas, Characterizing the inf-sup condition on product spaces, Numer. Math. 109 (2) (2008) 209–231. [22] R.H.W. Hoppe, P. Porta, Y. Vassilevski, Computational issues related to iterative coupling of subsurface and channel flows, Calcolo 44 (2007) 1–20. [23] O. Iliev, V. Laptev, On numerical simulation of flow through oil filters, Comput. Visual. Sci. 6 (2004) 139–146. [24] W. Jäger, A. Mikelic´ , On the interface boundary condition of Beavers, Joseph, and Saffman, SIAM J. Appl. Math. 60 (2000) 1111–1127. [25] M. Johnson, What controls aqueous humor outflow resistance? Exp. Eye Res. 82 (2006) 545–557. [26] H. Kolb, Gross anatomy of the eye, . [27] W. Layton, C. Manica, M. Neda, M.A. Olshanskii, L. Rebholz, On the accuracy of the rotation form in simulations of the Navier–Stokes equations, J. Comput. Phys. 228 (5) (2009) 3433–3447. [28] W.J. Layton, F. Schieweck, I. Yotov, Coupling fluid flow with porous media flow, SIAM J. Numer. Anal. 40 (2003) 2195–2218. [29] Y.-J. Lee, H. Li, On stability accuracy and fast solvers for finite element approximations of the axisymmetric Stokes problem by Hood–Taylor elements, SIAM J. Numer. Anal. 49 (2) (2011) 668–691. [30] H. Li, Finite element analysis for the axisymmetric Laplace operator on polygonal domains, J. Comput. Appl. Math. 235 (17) (2011) 5155–5176. [31] M. Mu, J. Xu, A two-grid method of a mixed Stokes–Darcy model for coupling fluid flow with porous media flow, SIAM J. Numer. Anal. 45 (5) (2007) 1801– 1813. [32] F.W. Newell, Ophthalmology: Principles and Concepts, volume 15 of Texts in Applied Mathematics, Mosby, St. Louis, 1996. eighth ed.. [33] M. Olshanskii, G. Lube, T. Heister, J. Löwe. Grad-div stabilization and subgrid pressure models for the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 198 (2009) 3975–3988. [34] M.A. Olshanskii, A low order Galerkin finite element method for the Navier– Stokes equations of steady incompressible flow: A stabilization issue and iterative methods, Comput. Methods Appl. Mech. Engrg. 191 (2002) 5515– 5536. [35] M.A. Olshanskii, A. Reusken, Grad-Div stabilization for the Stokes equations, Math. Comput. 73 (2004) 1699–1718. [36] International Commission on Radiological Protection, Task Group on Reference Man, Report of the Task Group on Reference Man, ICRP No. 23. Pergamon Press, Oxford, New York, 1975. [37] P.-A. Raviart, J.M. Thomas, Dual finite element models for second order elliptic problems, Energy Methods in Finite Element Analysis, Wiley, Chichester, 1979, pp. 175–191. [38] B. Riviére, I. Yotov, Locally conservative coupling of Stokes and Darcy flows, SIAM J. Numer. Anal. 42 (2005) 1959–1977. [39] P. Saffman, On the boundary condition at the surface of a porous media, Stud. Appl. Math. 50 (1971) 93–101.