Hybridizable discontinuous Galerkin method (HDG) for Stokes interface flow

Journal of Computational Physics 247 (2013) 262–278

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

Hybridizable discontinuous Galerkin method (HDG) for Stokes interface flow Bo Wang a, B.C. Khoo b,c,⇑ a

Singapore-MIT Alliance, 4 Engineering Drive 3, National University of Singapore, Singapore 117576, Singapore Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore c Temasek Laboratories, 5A Engineering Drive 1, Singapore 117411, Singapore b

a r t i c l e

i n f o

Article history: Received 4 August 2012 Received in revised form 11 March 2013 Accepted 26 March 2013 Available online 12 April 2013 Keywords: Stokes interface problem Hybridizable discontinuous Galerkin method Curvilinear element Postprocessing

a b s t r a c t In this paper, we present a hybridizable discontinuous Galerkin (HDG) method for solving the Stokes interface problems with discontinuous viscosity and variable surface tension. The jump condition of the stress tensor across the interface is naturally incorporated into the HDG formulation through a constraint on the numerical flux. The most important feature of HDG method compared to other DG methods is that it reduces the number of globally coupled unknowns significantly when high order approximate polynomials are used. For problems with polygonal interfaces, it provides optimal convergence rates of order k þ 1 in L2 -norm for the velocity, pressure and as well as the gradient of velocity. Furthermore, a new approximate velocity can be obtained by an element-by-element postprocessing which converges with order k þ 2 in the L2 -norm. For Stokes interface problems with curved interfaces, we use general curvilinear element to ensure the optimal convergence rates. An error estimate is given for the approximation of the interface. It indicates that curvilinear elements of degree 2k þ 1 should be used for optimal convergence rate of order k þ 1. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction Stokes interface problems arise from many applications such as fluid structure interaction (FSI) problems [12,17,30,31], two-phase incompressible flows [3,10,18,19] and so on. In FSI problems, an interface force is exerted on the boundary of the structure which results in a jump in the flux across the boundary. In two-phase incompressible flow, fluids with different density and viscosity are involved. The density and viscosity are then discontinuous across the interface between the two different fluids. Therefore, a jump condition of the flux across the interface needs to be taken into account in dealing with two-phase incompressible flows. In the implementation for low Reynolds number flow regime, a Stokes interface problem needs to be solved at each time step. One popular numerical method for solving FSI problem is the immersed boundary method, see [27–29]. The main idea of this method is to model the effect of the interface force by an equivalent singular force on the interface and then redistribute this singular force to the vicinity of the interface using a so-called regularized delta function. This idea also has been coupled with level set method for solving two phase incompressible flows, see [3,20,25]. The most important feature of this strategy is that the mesh do not need to fit the interface. Therefore, there is no need to re-mesh the domain when the interface

⇑ Corresponding author at: Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore. Tel.: +65 68742889. E-mail addresses: [email protected] (B. Wang), [email protected] (B.C. Khoo). 0021-9991/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.03.064

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

263

moved. However, the redistribution of the singular force results in first order accuracy even when high order numerical schemes are used in the discretization. Another important numerical method for solving interface problems is the immersed interface method (IIM). The IIM has been widely used in interfacial elliptic problems, Stokes problems and Navier–Stokes problems [15–19,30,31], just to name a few. This method incorporate the known jumps in the solutions or their derivatives into the finite difference scheme to obtain a modified scheme. It is second order accurate at all points in the uniform grid. Some high order IIM are also available for special interface problems, see [14]. In the construction of the formulation using IIM, there is the need of jump conditions for the solutions and their derivatives. When dealing with FSI or two phase flow problems one usually has only the jump condition for the velocity and normal stress pn þ mru. Therefore, in order to use IIM, one would require to derive the necessary jump conditions for the velocity, pressure and their derivatives based on the model problem and known jump conditions. This procedure can and usually is (fairly) complex especially for the 3D problems. It is worthwhile to mention that a theory on uniform regularity for the Stokes interface problem has been investigated in [26]. Uniform regularity results of the velocity and pressure with respect to the jump in the viscosity has been proved. Furthermore the authors use a LBB-stable finite element method to solve the Stokes interface problem on a conforming mesh. A robust preconditioner for the Schur complement is given and shown to be very efficient theoretically and also numerically. Recently, the hybridizable discontinuous Galerkin (HDG) method is developed and used to solve the elliptic, Stokes, Navier–Stokes, and Maxwell equations [6,5,23,24,22,21]. The HDG method is very attractive due to the reduction in the number of globally coupled unknowns especially for the high order case. The hybridization of the DG methods for Stokes equations was initially introduced in [2] to get a divergence-free velocity without actually constructing a divergence-free finite element space for the approximation of the velocity. Based on the hybridized technique, the first HDG methods for the Stokes equations were introduced in [5]. The hybridization in four completely different ways are discussed in the said paper. In [23], a HDG method based on velocity–pressure-gradient formulation is introduced. The approach maintains the divergence free property of velocity. Numerical solutions including velocity, pressure and gradient of velocity exhibited an optimal order of k þ 1 in the L2 -norm, when polynomials of degree k are used as the basis of finite element spaces. Further, a new approximate velocity which converges with order k þ 2 in the L2 -norm can be obtained by an element-by-element postprocessing. In this paper, we investigate the HDG method based on the velocity–pressure-gradient formulation for solving the Stokes interface problems in a 2-D polygonal domain. The jump condition of the stress tensor is naturally incorporated into the HDG formulation through a constraint on the corresponding numerical flux. It is therefore not necessary to derive any additional jump conditions for the velocity, pressure and their derivatives as for the IIM. If the interface is polygonal, a linear conforming mesh is used. We show that the velocity, pressure and as well as the gradient of velocity have uniform optimal convergence rate of order k þ 1 in the L2 -norm with respect to the jump in the viscosity. However, linear conforming mesh can not maintain optimal convergence rate when curved interface is involved. In fact, the significant influence of the geometry representation of curved interface or boundary on the accuracy of numerical methods has already been noticed and investigated, see [11,9] and their references. As such, in order to obtain optimal convergence rate we use the curvilinear element in the discretization. After the approximation of the interface we obtain an approximate problem. An error estimate is given for this approximation based on the so-called continuum surface force (CSF) model. It indicates that curvilinear element of degree 2k þ 1 should be used for optimal convergence rate of order k þ 1. Upon the optimal convergence rate for the velocity and the gradient of velocity in L2 -norm, a new approximate velocity with convergence rate of order k þ 2 can be obtained by using element-by-element postprocess. The paper is organized as follows. In Section 2, we introduce the model problem and its variational formulation. In Section 3, the curvilinear element is introduced. An error estimate for the approximation of the interface is given. The HDG formulation is then presented in Section 4. We show how the jump conditions are naturally incorporated into the formulation. Two implementation strategies are discussed in detail. After establishing the HDG scheme, we present a L2 -error estimate. Using the augmented Lagrangian implementation strategy we present some numerical results in Section 5. Finally, concluding remarks are given in Section 6. 2. Model problem In this paper we consider the following interface Stokes problem

r  ðmðxÞruÞ þ rp ¼ f

in Xk ;

r  u ¼ 0 in Xk k ¼ 1; 2;

ð1aÞ ð1bÞ

sutC ¼ 0;

srðu; pÞtC n ¼ g;

ð1cÞ

uj@ X ¼ ub

on @ X;

ð1dÞ

with a piecewise constant viscosity

m¼



m1 in X1 ; m2 in X2 ;

ð2Þ

264

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

on a bounded connected polygonal domain X. And X1 and X2 are Lipschitz subdomains of X such that X1 \ X2 ¼ ; and 1 [ X  2 . The interface between the subdomains is denoted by C ¼ @ X1 \ @ X2 . See Fig. 1 for the configuration of the do ¼X X main. Further, rðu; pÞ is the stress tensor which is given by rðu; pÞ ¼ pI þ mru and n is a unit normal vector to C pointing from X1 to X2 . The jump notation is given by sv tC ¼ ðv jX1  v jX2 ÞjC . It is well known that the boundary data ub must satisfy the compatibility condition

Z

ub  n ¼ 0;

ð3Þ

@X

due to the incompressibility constraint. Here, for simplicity of analysis, we shall assume homogeneous boundary conditions, i.e, ub ¼ 0. However, non-homogenous boundary conditions will be considered in our numerical examples. To determine p R uniquely we may impose some additional condition, here we choose X p ¼ 0. k We shall use the following spaces. Denoted by H ðDÞ the standard Sobolev space in domain D and its norm denoted by P k  kHk ðDÞ . We set the corresponding vector space Hk ðDÞ :¼ ½Hk ðDÞn and kv kHk ðDÞ :¼ ni¼1 kv i kHk ðDÞ ; and tensor space P n nn Hk ðDÞ :¼ ½Hk ðDÞ and kZkHk ðDÞ :¼ i;j¼1 kZ ij kHk ðDÞ . Let L2 ðDÞ be the space of square integrable functions on D. Further denote their corresponding vector and tensor spaces by L2 ðDÞ ¼ ½L2 ðDÞd , and L2 ðDÞ ¼ ½L2 ðDÞdd . Define spaces

n o  Hk0 ðDÞ :¼ u 2 Hk ðDÞuj@D ¼ 0 ;

L20 ðDÞ :¼

Z    p 2 L2 ðDÞ pðxÞdx ¼ 0 : D

H10 ð

1

Set H ðXÞ be the dual space of

XÞ with the associated norm: jhf; v ij

f 2 H1 ðXÞ # kfkH1 ðXÞ :¼ sup

v2H10 ðXÞ kv kH1 ðXÞ

The two Stokes equations in (1a) and the coupling jump condition across the interface in (1c) can be reformulated into one Stokes equation in the whole domain in which the effect of the surface tension is expressed in terms of a localized force at the interface; cf. the so-called continuum surface force (CSF) model [1]. We define the variational formulation of the problem (1) as follows: find ðu; pÞ 2 H10 ðXÞ  L20 ðXÞ such that

ðmru; rv Þ  ðr  v ; pÞ ¼ ðf; v Þ þ f C ðv Þ for ðr  u; qÞ ¼ 0 for q 2

v 2 H10 ðXÞ;

L20 ð

XÞ;

ð4aÞ ð4bÞ

with

f C ðv Þ ¼

Z

srðu; pÞtn  v ds ¼ C

Z

g  v ds:

ð5Þ

C

Here, f C is a functional on H10 ðXÞ if g 2 ½L2 ðCÞ2 . According to the standard theory of Stokes equations, we have: Lemma 1. Assume that f 2 ½L2 ðXÞ2 and g 2 ½L2 ðCÞ2 ; m > 0 is piecewise constant, then the variational problem (4) has a unique solution and that the priori estimate 1

ðkuk2H1 ðXÞ þ kpk2L2 ðXÞ Þ2 6 CðkfkH1 ðXÞ þ kf C kH1 ðXÞ Þ holds with a generic constant C depends on m.

Fig. 1. A diagram of the geometry of the computational domain.

ð6Þ

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

265

3. Curvilinear elements for curved interfaces In this section, we introduce the meshes we used for problems with polygonal interface and curved interface, respectively. Curvilinear elements are used for problems with curved interface to ensure the optimal convergence rates. We shall give an error estimate for the approximated problem resulting from approximating the interface by curvilinear elements. Denote by T h a collection of disjoint regular elements K that partition X and set @T h :¼ f@K : K 2 T h g. For an element K of the collection T h ; F ¼ @K \ @ X is the boundary face if the d  1 Lebesgue measure of F is nonzero. We shall define the mesh T h as a conforming mesh w.r.t. the two subdomains X1 ; X2 in the following sense: ðiÞ

9T h  T h :

o [n ðiÞ  i; KjK 2 T h ¼ X

i ¼ 1; 2;

ð7Þ

where X1 and X2 are polyhedral subdomains. As such, a conforming mesh can be obtained when the interface is polygonal. We use this conforming mesh when the interface is polygonal. However, in most computational fluid dynamic problems it is (more) realistic to assume that the interface C is a smooth curve. In this paper, our interest is in the Stokes interface problems in which the interface C is curved. In solving such problems, the geometric representation of the interface C has a significant influence on the accuracy of numerical method [11]. We shall use curvilinear elements to improve the accuracy of the approximation of the interface so that optimal convergence rates can be obtained. The curvilinear elements are based on a linear conforming mesh. In a linear conforming mesh, the interface intersects with the mesh only on the nodal points. Hence the interface C is divided into N parts fCi gNi¼1 by some nodal points in linear conforming mesh, see Fig. 2 (the interface is divided into 8 parts). In order to obtain a more accurate geometric representation of the interface we need to curve the edges whose end points are located on the interface according to the exact description of the curved interface. This can be done by using high order polynomials to approximate each Ci . As such, for each Ci we obtain a sufficiently accurate approximation Cih , see Fig. 3 (polynomials of degree 2 are used to approximate a Ci ). We replace the linear edge by the curvilinear edge Cih to obtain a curvilinear element. By this way, the two subdomains X1 and X2 are approximated by Xh;1 and Xh;2 , where

X ¼ Xh;1 [ Xh;2 :

Fig. 2. Linear conforming mesh (left) and orthogonal projection on C (right).

Fig. 3. A diagram of the curvilinear element.

266

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

It is clear that the curvilinear element approximate the curved interface much better than the linear element. Furthermore, we have the flexibility to use high order polynomials to obtain the approximation with desired accuracy. We need to introduce some notations in Fig. 3 which will be used in the error estimate. We see that there are some small domains bounded by i2 Ci and Cih . These small domains are denoted as Aijh , e.g. two small domains (shadowed area) Ai1 h and Ah in Fig. 3. Set

Aih ¼

MðiÞ [

Aijh ;

Ah ¼

N [ Aih ;

ð8Þ

i¼1

j¼1

where MðiÞ is the number of small domains attaching Ci and Cih . N S Denoted by Ch ¼ Cih the approximated interface of C. Since the jump of the stress tensor srðu; pÞtC is only defined on i¼1

the interface C, we need to extend it to Ch . For simplicity, we denote srðu; pÞtC by Jr . In this paper, we extend Jr to Ah in the following way:

Jr ðx; yÞ ¼ Jr ðx ; y Þ for ðx; yÞ 2 Aih 

i ¼ 1; 2; . . . ; N;

ð9Þ i





where ðx ; y Þ is the orthogonal projection of the point ðx; yÞ on the interface C . In Fig. 2, P is the orthogonal projection of P on C, hence Jr ðPÞ ¼ Jr ðP Þ. Here, we have to make sure that all points in Aih have unique orthogonal projection on Ci . This constraint is reasonable. For interface with large curvature we can use fine mesh near the interface. After the approximation of the interface, we obtain a new Stokes interface problem:

~ Þ þ rp ~ ¼ f in Xh;k ;  r  ðmðxÞru ~ ¼ 0 in Xh;k k ¼ 1; 2; ru ~ tC ¼ 0; srðu ~; p ~ÞtC n ¼ Jr n; su

ð10bÞ

~ j@ X ¼ ub u

ð10dÞ

h

ð10aÞ ð10cÞ

h

on @ X;

where, n is a unit normal vector to Ch pointing from Xh;1 to Xh;2 . The variational form of (10) is: Determine ~; p ~Þ 2 H10 ðXÞ  L20 ðXÞ such that ðu

~ ; rv Þ  ðr  v ; p ~Þ ¼ ðf; v Þ þ f Ch ðv Þ for ðmru ~ ; qÞ ¼ 0 for q 2 ðr  u

v 2 H10 ðXÞ;

L20 ð

XÞ;

ð11aÞ ð11bÞ

with continuum surface force

f Ch ¼

Z Ch

Jr n  v ds:

ð12Þ

This is the problem which the HDG method is used to solve in this paper. We need to point out that the viscosity m is not changed due to the approximation of the interface. Only the surface tension is replaced by an approximate singular force on Ch . Therefore, we could give the following estimate between the solutions of problem (4) and (11). ~; p ~Þ be the solution of (4) and (11), respectively. Then Theorem 2. Let ðu; pÞ and ðu 1

~ k2H1 ðXÞ þ kp  p ~k2L2 ðXÞ Þ2 6 CkJr kH1 ðA Þ : ðku  u h

ð13Þ

~u ¼ u  u ~ and ~ ~, we have Proof 1. Denoted by e ep ¼ p  p

~u ; rv Þ  ðr  v ; ~ep Þ ¼ f C ðv Þ  f Ch ðv Þ for ðmre ~u ; qÞ ¼ 0 for q 2 ðr  e

v 2 H10 ðXÞ;

L20 ð

XÞ:

ð14aÞ ð14bÞ

Lemma 1 indicates that we have the estimate 1

~u k2H1 ðXÞ þ k~ep k2L2 ðXÞ Þ2 6 Cðkf C  f Ch kH1 ðXÞ Þ: ðke

ð15Þ

Since we use curvilinear elements to approximate the interface, C and Ch must form a series of sub-domains, see Fig. 3. Divergence theorem shows that

     Z Z MðiÞ Z  X  X   N MðiÞ Z N X   X     T T ~ ds ¼  jf C ðv Þ  f Ch ðv Þj ¼  Jr n  v ds  Jr n  v ds ¼  Jr v  n r  ðJr v Þdx ij ij     C   Ch i¼1 j¼1 @Ah i¼1 j¼1 Ah    X N Z   ¼ ðr  Jr  v þ Jr : rv Þdx:   i¼1 Ai h

~ is outward normal of Aijh on its boundary. Then we have Here, n

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

jf C ðv Þ  f Ch ðv Þj 6

Z

267

ðjr  Jr kv j þ jJr krv jÞdx

Ah

Z 6 Ah

!12 Z

jr  Jr j2 dx

!12 jv j2 dx

Z

þ

Ah

Ah

jJr j2 dx

!12 Z

!12 jrv j2 dx

Ah

6 ðkr  Jr kL2 ðAh Þ þ kJr kL2 ðAh Þ Þkv kH1 ðXÞ ;

ð16Þ

by using Schwartz inequality. This implies

kf Ch  f C kH1 ðXÞ 6 kJr kH1 ðAh Þ :

ð17Þ



Remark 1. We would like to point out that this estimate is available for other types of approximation of the interface. We note that the regularity of function Jr ðx; yÞ in Ah depends on the function srðu; pÞtC and the regularity of interface as well. Assume that x ¼ XðsÞ; y ¼ YðsÞ is the arc-length parametrization of the interface C. According to the definition of the extension (9) we have

 @Jr  ¼ 0; @nP P

 @Jr  dsrðu; pÞtC ¼ @ gP P ds

ð14rÞ

where nP and gP are the normal and tangential vector at the point P . We can solve rJr from (14r) giving rise to

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2  2 @XðP Þ @YðP Þ rJr ðPÞ ¼ þ @s @s

@YðP  Þ @s @XðP  Þ @s

 @XðP @s



Þ

!1

@YðP Þ @s

0 dsrðu;pÞtC ds

! :

ð14sÞ

Theorem 3. Assume that the interface is piecewise in C 1 , and

  dsrðu; pÞtC    6 CJ;   ds

jsrðu; pÞtC j 6 C J ;

then we have 1

~ k2H1 ðXÞ þ kp  p ~k2L2 ðXÞ Þ2 6 CC J ðku  u

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi areaðAh Þ;

ð14tÞ

where C is a generic constant.  If we use polynomials of degree at most k in the construction of the curvilinear elements, we will have the estimate, k þ1 areaðAh Þ 6 Ch . Hence, for ensuring optimal convergence rate of order k þ 1 we have to use polynomials of degree at most  k ¼ 2k þ 1 to approximate the interface. This type of elements is also named as super-parametric elements [13]. 4. HDG formulation 4.1. Some notations For the description of the HDG method we shall introduce some notations as used in [23,22]. For two elements K þ and K  of the mesh collection T h ; F ¼ @K þ \ @K  is the interior face between K þ and K  if the d  1 Lebesgue measure of F is nonzero. Denote by E oh and E @h the set of interior and boundary faces, respectively. We set E h ¼ E oh [ E @h . Let nþ and n be the outward unit normal vectors on two neighboring elements K þ and K  , respectively. We use ðG ; v  ; q Þ to denote the trace of ðG; v ; qÞ on F from the interior of K  , where G; v and q are second-order tensorial, vectorial and scalar functions, respectively. For interior face F 2 E oh , we set

sGnt ¼ Gþ nþ þ G n ;

sv nt ¼ v þ nþ þ v  n ;

sqnt ¼ qþ nþ þ q n

and for boundary face F 2 E @h on which G; v and q are single-valued, we set

sGnt ¼ Gn;

sv nt ¼ v n;

sqnt ¼ qn:

Here is either  or which denote the usual dot product and tensor product, respectively. Next we introduce the finite element spaces. Let P k ðDÞ denote the space of polynomial of degree at most k on a domain D. Further denote their corresponding vector and tensor spaces by Pk ðDÞ ¼ ½P k ðDÞd ; Pk ðDÞ ¼ ½P k ðDÞdd . We introduce the following discontinuous finite element approximation spaces for the gradient, velocity, and pressure:

Gh ¼ fG 2 L2 ðT h Þ : GjK 2 Pk ðKÞ; 8K 2 T h g; Vh ¼ fv 2 L2 ðT h Þ : v jK 2 Pk ðKÞ; 8K 2 T h g; Ph ¼ fq 2 L2 ðT h Þ : qjK 2 P k ðKÞ; 8K 2 T h g:

268

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

A finite element approximation space Mh for trace of the velocity is defined as

Mh ¼ fl 2 L2 ðE h Þ : ljF 2 Pk ðFÞ; 8F 2 E h g and denoted by Mh ðub Þ ¼ fl 2 Mh : l ¼ P@ ub on @ Xg, where P@ is the L2 ð@ XÞ projection into the space flj@ X ; 8l 2 Mh g. The  h and defined as mean of the approximate pressure belongs to a space denoted by W

 h ¼ fr 2 L2 ð@T h Þ : r 2 P 0 ð@KÞ; 8K 2 T h g: W This space consists of functions that are constant on each @K for all element K 2 T h . The mean of a function q 2 L2 ð@T h Þ on @K R 1  h. j@K ¼ j@Kj  ¼ q for all q 2 W is given by q q. Obviously, we have q @K Finally, we shall define various inner products for our finite element spaces. Define

X

ðu; v ÞT h :¼

ðu; v ÞK ;

K2T h

where ðu; v ÞD denotes the integral of uv over the domain D  Rn . Further, define

ðu; v ÞT h :¼

n X ðui ; v i ÞT h

and ðN; ZÞT h :¼

i¼1

n X ðNij ; Z ij ÞT h i;j¼1

and

X

hu; v i@T h :¼

hu; v i@K

and hu; v i@T h :¼

K2T h

n X hui ; v i i@T h ; i¼1

where hu; v iD denotes the integral of uv over the domain D  Rn1 . 4.2. The formulation In order to devise the HDG formulation, we reformulate (10) in velocity–pressure-gradient form as follows:

L  ru ¼ 0;

in X n Ch ;

 r  rðL; pÞ ¼ f; in X n Ch ; r  u ¼ 0; in X; sutCh ¼ 0; u ¼ ub ;

srðL; pÞtCh n ¼ Jr n;

on @ X;

ð21aÞ ð21bÞ ð21cÞ ð21dÞ ð21eÞ

~ and p ~ in (10) by u and p in (21) and the rest of the paper for where rðL; pÞ ¼ pI þ mL. Note that we have denoted u simplicity. ~ ¼ Jr njC and introduce the HDG formulation for system (21): Find an approximation ðLh ; uh ; ph ; u ^ h Þ of Next, we denote g h ^ jE Þ in the space Gh  Vh  P h  Mh such that ðLjX ; ujX ; pjX ; u h

^ h ; Gni@T ¼ 0; ðLh ; GÞT h þ ðuh ; r  GÞT h  hu h

ð22aÞ

^ h Þn; v i@T ¼ ðf; v ÞT ; ^ ðLh ; uh ; ph ; u ðrðLh ; ph Þ; rv ÞT h þ hr h h

ð22bÞ

^ h  n; qi@T ¼ 0;  ðuh ; rqÞT h þ hu h

ð22cÞ

^ h Þn; li@T n@ X ¼ hg ~ ; li C ^ ðLh ; uh ; ph ; u hr h h

ð22dÞ

^ h ; li@ X ¼ hub ; li@ X ; hu

ð22eÞ

ðph ; 1ÞT h ¼ 0

ð22fÞ

for all ðG; v ; q; lÞ 2 Gh  Vh  P h  Mh . The numerical flux has the form

r^ ðLh ; uh ; ph ; u^ h Þ ¼ ph I þ mLh  Sðuh  u^ h Þ n:

ð23Þ

Here, S is a second order tensor consisting of stabilization parameter. It plays a key role on the stability and convergence of HDG method. Some general requirements on the stabilized tensor S are proposed in [4]. In this paper, we use a simple special choices of S, i.e.,

( S¼

msI

on F H K;

0

on @K n F H K;

ð24Þ

where F H K is an arbitrary face of K. Obviously, this special choice must satisfy all general requirements as stated in [4] with

kmin ¼ Kmax ¼ ms: K K

269

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

Except for the forth Eq. (22d), formulation (22) is the same as that of HDG method for conventional Stokes problem presented in [4]. In fact, the change in Eq. (22d) indicates that the jump condition of the flux are naturally imposed in the weak sense, i.e.

Z

^ h Þnt  l ds ¼ ^ ðLh ; uh ; ph ; u sr

F

Z

~  l ds; g

8F 2 C h :

ð25Þ

F

More clearly, we have

Z F

ðsrðLh ; ph Þnt  ssuh tÞ  l ds ¼

Z

~  l ds; g

8F 2 C h

ð26Þ

F

^ h is a single valued function defined on E h . for u The first three equation in (22) define the following local solver. The local solver Lh maps any given data  to the function ðLh ; uh ; p Þ 2 Gh  Vh  P h satisfying  2 L2 ðXÞ  Vh  W ðf; g; wÞ h

ðLh ; GÞK þ ðuh ; r  GÞK ¼ hg; Gni@K ;

ð27aÞ

ðrðLh ; ph Þ; rv ÞK þ hrðLh ; ph Þn þ Suh ; v i@K ¼ ðf; v ÞK þ hSg; v i@K ;

ð27bÞ

 ðuh ; rqÞK ¼ hg  n; qi@K ;

ð27cÞ

 h ¼ w: p

ð27dÞ

We shall denote the following:

ðLfh ; ufh ; pfh Þ :¼ Lh ðf; 0; 0Þ;

ð28aÞ

ðLgh ; ugh ; pgh Þ :¼ Lh ð0; g; 0Þ;

ð28bÞ

    ðLwh ; uwh ; pwh Þ :¼ Lh ð0; 0; wÞ:

ð28cÞ









 in the local solver Remark 2. It is worthwhile to point out that Lwh ¼ 0 and uwh ¼ 0. In fact, setting G ¼ mLwh ; v ¼ uwh and q ¼ w we have 







mðLwh ; Lwh ÞK þ hSuwh ; uwh i@K ¼ 0; 8K 2 T h ; which indicates the fact.  is eliminated from the third Note that there is a small difference from the local solver defined in [23]; here the quantity q Eq. (27c) of the local solver. However, this change will not affect the formulation. In fact, by using the same technique as in [23] we still have a similar theorem. ^ h Þ be the solution of (22). We have that Theorem 4. Let ðLh ; uh ; ph ; u

Lh ¼ Lfh þ Lkh uh ¼ ufh þ ukh ph ¼ pfh þ pkh þ pqh ^h ¼ k u

ð29Þ

h ¼ q p  h Þ satisfies where ðk; qÞ 2 ðMh ðub Þ; W

ah ðk; lÞ þ bh ðq; lÞ ¼ lh ðlÞ; h  kÞ ¼ 0; 8w  2W bh ðw;

8l 2 Mh ð0Þ;

ð30Þ

and

ðpfh þ pkh þ pqh ; 1ÞT h ¼ 0:

ð31Þ

Here, the forms are given by l

l

ah ðg; lÞ ¼ ðmLgh ; Lh ÞT h þ hSðugh  gÞ; ðuh  lÞi@T h ;  lÞ ¼ hw  h ; l  ni ; bh ðw; @T h

ð32Þ

l

~; liC lh ðlÞ ¼ ðf; uh ÞT h þ hg h  h.  2W for all g; l 2 Mh and w Here, the linear form lh ðlÞ in Theorem 4 is different from that introduced in [23] due to the weakly imposition of the jump condition on the flux.

270

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

The weak formulation (30) in Theorem 4 give rise to a linear system of the form

A BT B

!

0

K !

 ¼

  F ; 0

ð33Þ

where K and ! represent the vectors of freedom of k and q, respectively. Once we get the approximate trace of velocity k and  ¼ q. the mean value of pressure q, the HDG solution can be obtained by solving (27) element-by-element with g ¼ k and w There is another efficient way to implement the HDG method i.e. the augmented Lagrangian method. Giving an initial pressure p0 2 L20 ðXÞ, construct an evolution problem

@pðt  Þ þ r  uðt  Þ ¼ 0; in X  ð0; 1Þ; @t  pðt  ¼ 0Þ ¼ p0 ; on @ X; 

ð34aÞ ð34bÞ



where uðt Þ is a function of pðt Þ and defined as the solution of

Lðt  Þ  ruðt  Þ ¼ 0;

in X  ð0; 1Þ;

ð35aÞ

 r  rðL; pÞ ¼ f;

ð35bÞ

sutCh ¼ 0;

ð35cÞ



in X  ð0; 1Þ; ~; srðL; pÞntCh ¼ g

uðt Þ ¼ ub ;

on @ X  ð0; 1Þ:

ð35dÞ 



The system (34) and (35) is the time-dependent version of the original problem (21). It can be proved that Lðt Þ; uðt Þ and R pðt Þ, the solutions of (34) and (35), converge exponentially to the solutions of (21) with restriction X p ¼ 0 as the pseudo time t goes to infinity. Thus the HDG scheme based on the augmented system (34) and (35a)–(35d) can be defined as follows: given an initial guess

ðp0h ; qÞT h ¼ ðp0 ; qÞT h

ð36Þ

and the iteration procedure is given by

^ nh ; Gni@T ¼ 0; ðLnh ; GÞT h þ ðunh ; r  GÞT h  hu h ðr

ðLnh ; pnh Þ;

rv ÞT h

^ nh Þn; v i@T ¼ ðf; v ÞT ; ^ ðLnh ; unh ; pnh ; u þ hr h h

ð37aÞ ð37bÞ

1 1 ^ nh  n; qi@T ¼  ðpn1 ðpn ; qÞ  ðunh ; rqÞT h þ hu ; qÞT h ; h Dt  h T h Dt h n ^ nh Þn; li@T ¼ hg ~; liC ^ ðLh ; unh ; pnh ; u hr h h

ð37dÞ

^ h ; li@ X ¼ hub ; li@ X ; hu

ð37eÞ

ð37cÞ

for all ðG; v ; q; lÞ 2 Gh  Vh  P h  Mh . Similarly, we can see that (37a)–(37c) define a local solver. We denote it by Lh mapping any given data ðf; g; nÞ 2 L2 ðXÞ  Vh  Ph to the function ðLh ; uh ; ph Þ 2 Gh  Vh  P h satisfying

ðLh ; GÞK þ ðuh ; r  GÞK ¼ hg; Gni@K ;

ð38aÞ

ðrðLh ; ph Þ; rv ÞK þ hrðLh ; ph Þn; v i@K þ hSuh ; v i@K ¼ ðf; v ÞK þ hSg; v i@K ;

ð38bÞ

1 1 ðp ; qÞ  ðuh ; rqÞK ¼ hg  n; qi@K þ  ðn; qÞK : Dt  h T h Dt

ð38cÞ

Denoted by

ðLfh ; ufh ; pfh Þ :¼ Lh ðf; 0; 0Þ;

ð39aÞ

ðLgh ; ugh ; pgh Þ :¼ Lh ð0; g; 0Þ;

ð39bÞ

ðLnh ; unh ; pnh Þ :¼ Lh ð0; 0; nÞ;

ð39cÞ

we then have the following results, see [23] for details. ^ nh Þ be the solution of the system (37), then we have Lemma 5. If ðLnh ; unh ; pnh ; u pn1

Lnh ¼ Lfh þ Lkh þ Lhh

pn1

unh ¼ ufh þ ukh þ uhh

pn1

pnh ¼ pfh þ pkh þ þphh ^h ¼ k u

ð40Þ

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

271

where k 2 Mh ðub Þ satisfies

ch ðk; lÞ ¼ fh ðl; pn1 h Þ;

8l 2 Mh ð0Þ:

ð41Þ

Here the forms are given by l

l

ch ðg; lÞ ¼ ðmLgh ; Lh ÞT h þ hSðugh  gÞ; ðuh  lÞi@T h þ l

fh ðlÞ ¼ ðf; uh ÞT h 

1 l ðpg ; p Þ ; Dt  h h T h

1 l ~; liC ðpn1 ; ph ÞT h þ hg h Dt  h

for all g; l 2 Mh . Therefore, for each iteration we need to solve the following linear system

CK ¼ F;

ð42Þ

derived from the formulation (41). Here, K represents the vector of freedom of the numerical trace of velocity. Obviously, matrix C is symmetric positive definite if we impose the boundary condition appropriately. Hence, conjugated gradient (CG) method or preconditioned conjugated gradient (PCG) method can be used to solve (42). Once we obtain k, the HDG solution ðLnh ; unh ; pnh Þ can be obtained by solving (38) element-by-element with g ¼ k; n ¼ pn1 h . We stop the iteration when the relative error of pressure is less than a given tolerance tol , i.e.,

kpNh  pN1 kT h h kpNh kT h

< tol :

ð43Þ

Here, we briefly comment/summarize on the configuration of the two implementation strategies of the HDG method. In the first implementation strategy we need to solve a saddle point problem (33). An CG-Uzawa method may be employed. In the second implementation strategy, an iteration is needed, however, only the freedoms of the trace of velocity are involved in the global system (42) which permits it to have much less unknowns than that in the system (33). Furthermore, the coefficient matrix C in the global system (42) can be symmetric positive definite if the boundary condition is imposed appropriately. Thus, conjugated gradient (CG) method or preconditioned conjugated gradient (PCG) method can be used to solve (42). In our numerical experiments the second implementation strategy will be employed. 4.3. Local postprocessing A better velocity approximation uh can be obtained by a local postprocessing [23]. Here we show that this local postprocessing is also available for Stokes interface problem. By using curvilinear elements we approximate the curved interface very accurately. This helps us to obtain optimal convergence rate for both uh and Lh . Then we find uh where uh jK 2 P kþ1 ðKÞ such that

ðruh ; rv ÞK ¼ ðLh ; rv ÞK ; 8v 2 Pkþ1 ðKÞ

ð44Þ

ðuh ; 1ÞK ¼ ðuh ; 1ÞK

ð45Þ

for all K 2 T h . Apparently, the construction of uh is done element by element therefore it’s very efficient. The obtained new velocity approximation uh which has a convergence rate of order k þ 2 in L2 -norm as shown in our numerical results. 4.4. Error estimates If the interface is polygonal and a conforming mesh is used we have the following error estimates: Theorem 6. Suppose that s is positive, L; u and p is the exact solution of the problem (1), Lh ; uh and ph is the solution of (22). Then

kL  Lh kK 6 Ch

  1  L  pI  

kþ1 

kþ1

kp  ph kK 6 CC p h ku  uh kK 6 Ch

kþ1

m

; Hkþ1 ðKÞ

jmL  pIjHkþ1 ðKÞ ;

jujHkþ1 ðKÞ þ C

where

n pffiffiffiffiffiffio C p ¼ max 1; sh ; and C is a constant independent of m.

h

kþ1

s

    L  1 pI  m 

Hkþ1 ðKÞ

;

272

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

Denote the errors by EL ¼ L  Lh ; Eu ¼ u  uh , and Ep ¼ p  ph . It is obvious that the exact solution satisfies (22), hence we have the error equations ^

ðEL ; GÞT h þ ðEu ; r  GÞT h  hEu ; Gni@T h ¼ 0;

ð46aÞ ^

 ðr  ðmEL Þ; v ÞT h þ ðrEp I; v ÞT h þ hSðEu  Eu Þ; v i@T h ¼ 0;

ð46bÞ

^

 ðEu ; rqÞT h þ hEu  n; qi@T h ¼ 0; L

p

u

ð46cÞ

^ u

hmE n  E n  SðE  E Þ; li@T h n@ X ¼ 0;

ð46dÞ

^ u

hE ; li@T h \@ X ¼ 0;

ð46eÞ

ðph ; 1ÞT h ¼ 0;

ð46fÞ

for all ðG; v ; q; lÞ 2 Gh  Vh  Ph  Mh . Note that we have the same error equations as that for the regular Stokes equations. Therefore, the error estimate is just the same as that for the regular Stokes problems [4]. 5. Numerical examples In this section, some numerical examples are presented to validate the accuracy and fidelity of the approach. Example 1. In this example, we test a laminar flow in the domain X ¼ ½0; 2  ½0:5; 1:5, see Fig. 4. Two kinds of fluids whose viscosity are different flows in the subdomains X1 ¼ ½0; 2  ½0:5; 0:5, and X2 ¼ X n X1 , respectively. We keep the viscosity of the fluid in X1 fixed and vary the viscosity in X2 .That is, we set

m¼



m1 ¼ 1 in X1 ; m2 ¼  in X2

ð47Þ

and have the exact solution given as follows

u1 ¼ 1  ek sinð2pyÞ;

ð48Þ

u2 ¼ 0; 1 p ¼ e2kx ; 2

ð49Þ ð50Þ

where

8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 > < k1 ¼ Re1  Re1 þ 4p2 in X1 ; 2 4 k¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > : k ¼ Re2  Re22 þ 4p2 in X : 2 2 2 4

ð51Þ

Here, Re ¼ 1m is the Reynolds number. The pressure and the gradient of velocity are discontinuous across the interface y ¼ 0:5. The external force and boundary data can be evaluated from the exact solution. We compare the degree of the involved global freedom of HDG method with that of LDG method [7] in Table 1. In this table, N t and N F signifies the number of elements and edges respectively in the corresponding mesh. Dof is the degree of the global freedom. For the LDG method

Dof ¼ ½dimðPk ðKÞÞ þ dimðP k1 ðKÞÞ  Nt ;

ð52Þ

1.5

1 interface 0.5

0

−0.5 0

0.5

1

1.5

Fig. 4. Configuration of the computational domain.

2

273

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278 Table 1 1 Comparison of the number of global freedom for mesh with h ¼ 16 . Nt

NF

2048

k

3136

Dof

1 2 3 4

Percentage

LDG

HDG

14,336 30,720 53,248 81,920

14,592 20,864 27,136 17,728

101.79 67.92 50.96 21.64

while on the other hand,

Dof ¼ dimðP k ðFÞÞ  NF þ Nt

ð53Þ

for the HDG method. Table 1 shows that HDG method involves much less global freedom than LDG method when k increase. In order to show that the convergence rate of our HDG method is independent of the jump in viscosity, we test the cases with  ¼ 102 ; 104 . We choose s ¼ 1, the L2 -error and convergence rates are presented in Tables 2 and 3 for k ¼ 1; 2; 3, respectively. It is fairly clear that convergence rate of order k þ 1 is obtained for velocity, pressure and as well as the gradient of velocity. Further convergence rate of order k þ 2 is obtained for the postprocessing velocity uh . The augmented Lagrangian method is employed in our numerical tests. Table 4 shows the number of iterations N iter and corresponding CPU time when k ¼ 3 and 32  32 mesh is used. In this 32  32 mesh we have 2048 elements and 3136 edges. The table shows that N iter decrease as Dt  increases. However, the condition number of the linear system in each iteration increases as Dt increases. Therefore the pseudo time step size Dt is not the larger the better. Also we see that larger viscosity ratio mm12 increases the computation costs. We plot the numerical solution of the x-component of the velocity and the pressure with k ¼ 3 in Fig. 5. The numerical solution is calculated on 16  16 mesh and 64  64 mesh is used in the plot for showing the details in every element. For the purpose of showing the representative results, we compare the numerical solution on 16  16 mesh with the exact solution along the line x ¼ 1. Fig. 6 shows that the numerical velocity matches very well with the exact solution and the jump in pressure is well approximated by the HDG method. We also compare the numerical velocity uh with k ¼ 1 and its postprocessing velocity uh with the exact velocity along the line x ¼ 1 in Fig. 7. It is clear that the postprocessing velocity uh approximate the exact solution much better than uh .

Table 2 L2 -errors and convergence rates of uh ; ph ; Lh and uh with

 ¼ 102 .

k

h

ku  uh k

Order

kp  ph k

Order

kL  Lh k

Order

ku  uh k

Order

1.6889e1

–

1.3465e2

–

2.4905e1

–

1.4175e2

–

1

1 4 1 8 1 16

4.4760e2

1.9158

3.0895e3

2.1238

6.5589e2

1.9249

1.8092e3

2.9699

1.1396e2

1.9737

7.0695e4

2.1277

1.6738e2

1.9703

2.2839e4

2.9858

1 4 1 8 1 16

1.0888e2

–

4.0090e3

–

3.5840e2

1.4083e3

2.9507

6.4869e4

2.6276

4.4306e3

3.0160

7.3807e5

3.9704

1.7790e4

2.9848

8.7291e5

2.8936

5.4751e4

3.0165

4.6399e6

3.9915

1 4 1 8 1 16

2.1286e3

–

6.1025e4

–

4.1453e3

–

9.0549e5

–

1.3907e4

3.9360

4.9236e5

3.6316

2.6947e4

3.9418

2.9115e6

4.9589

8.8232e6

3.9784

3.3242e6

3.8886

1.7114e5

3.9783

9.2417e8

4.9775

2

3

Table 3 L2 -errors and convergence rates of uh ; ph ; Lh and uh with

1.1569e3

 ¼ 104 .

k

h

ku  uh k

Order

kp  ph k

Order

kL  Lh k

Order

ku  uh k

Order

2.4577e1

–

1.3433e2

–

3.6430e1

–

2.1024e2

–

1

1 4 1 8 1 16

6.5224e2

1.9138

3.0808e3

2.1244

9.5226e2

1.9357

2.6631e3

2.9809

1.6613e2

1.9731

7.0458e4

2.1285

2.4232e2

1.9744

3.3434e4

2.9937

1 4 1 8 1 16

1.6032e2

–

4.0093e3

–

5.2863e2

2.0707e3

2.9528

6.4870e4

2.6277

6.5313e3

3.0168

1.0861e4

3.9722

2.6144e4

2.9856

8.7292e5

2.8936

8.0690e4

3.0169

6.8240e6

3.9924

1 4 1 8 1 16

3.1120e3

–

6.1009e4

–

6.0770e3

–

1.3313e4

–

2.0198e4

3.9456

4.9229e5

3.6314

3.9269e4

3.9519

4.2671e6

4.9634

1.2779e5

3.9824

3.3232e6

3.8889

2.4867e5

3.9811

1.3527e7

4.9793

2

3

1.7046e3

274

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278 Table 4 The number of iterations and CPU time of the augmented Lagrangian method on 32  32 mesh with k ¼ 3; tol ¼ 1:0e  8.

 102

Dof

Dt 

N iter

CPU time (s)

25,088

8 16

18 12

8.7 6.18

32

9

4.92

8 16

18 12

25.37 17.51

32

9

17.72

25,088

104

ν =1, ν =1.0e−4

ν1=1, ν2=1.0e−4

2

2

0.6

1.5

0.4

1

0.2

p

p

1

0.5

0

0 −0.5

−0.2 −0.5

0 0.5 1 y

1.5

2

1

1.5

0.5

0

0

0.5 1

x

y

1.5

2

1

1.5

0.5

0

x

Fig. 5. Numerical solution of u and p on 16  16 mesh with k ¼ 3;  ¼ 104 .

In the next two examples, we shall consider a circular interface in a square domain. Since the interface is a smooth curve, we shall use curvilinear elements in the following numerical tests. Example 2. This example focuses on the normal force on an interface as in [8]. We shall consider a circular interface with radius 1 in a square domain X ¼ ½2; 22 . Let the interface force exerts only along the normal direction of the interface, and has the form

gðhÞ ¼ 2 sinð3hÞXðhÞ:

ð54Þ

Here XðhÞ ¼ ðcos h; sin hÞ is the parametrization of the interface. The exact pressure and velocity are given in polar coordinates:

( pðr; hÞ ¼ ( uðr; hÞ ¼ (

v ðr; hÞ ¼

r 3 sinð3hÞ; r < 1; r3 sinð3hÞ;

r P 1;

3 2 1 4 r sinð2hÞ þ 16 r sinð4hÞ  14 r 4 sinð2hÞ; 8 1 2 3 4 r sinð2hÞ  16 r sinð4hÞ þ 14 r 2 sinð4hÞ; 8 3 2 1 4 r cosð2hÞ  16 r cosð4hÞ  14 r 4 cosð2hÞ; 8 1 2 3 4 r cosð2hÞ þ 16 r cosð4hÞ  14 r2 cosð4hÞ; 8

r < 1; r P 1;

ð55Þ

r < 1; r P 1;

We note that p is discontinuous across the boundary, but ru; v and rv are continuous. The boundary data can be obtained from the exact solution on the boundary. And the jump condition is just the interface force given in (54) or one can compute by using the exact solution. The L2 -norm of the error of velocity, pressure and postprocessing velocity and convergence rate are presented in Table 5. We have optimal convergence rates of order k þ 1 for both velocity and pressure. The convergence rate is a little inferior to optimal order when k ¼ 3. In the next example, there is no such problem. Since we have a convergence rate of order k þ 2 for the postprocessing velocity, it serve to indicate that we have the optimal convergence rate for the gradient of velocity. We plot the numerical solution of the x-component of velocity 1 and the p with k ¼ 3 in Fig. 8. The numerical solution is calculated on the mesh with h ¼ 14 and mesh with h ¼ 16 is used in the plot for showing the details in every element. @p ; u; @r

275

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278 ν =1, ν =1.0e−4 1

2

ν1=1, ν2=1.0e−4

2

0.6

1.8 0.5

1.6 1.4

0.4 0.3

1

p

u

1.2

0.8

0.2

0.6

0.1

0.4 0

0.2 0 −0.5

0

0.5 x

1

−0.1 −0.5

1.5

0

0.5 x

1

1.5

Fig. 6. Numerical solution of u and p along the line x ¼ 1 on 16  16 mesh with k ¼ 3;  ¼ 104 (solid). Also plotted is the exact solution (star).

ν =1, ν =1.0e−4

ν1=1, ν2=1.0e−4

2.5

1

2

2

1.8 2

1.6 1.4

1.5

1

1

u*

u

1.2

0.8 0.5

0.6 0.4

0

0.2 −0.5 −0.5

0

0.5 x

1

0 −0.5

1.5

0

0.5 x

1

1.5

Fig. 7. Numerical solution of u and postprocessing velocity uh along the line x ¼ 1 on 16  16 mesh (solid) with k ¼ 1;  ¼ 104 . Also plotted is the exact solution (star).

Table 5 L2 -errors and convergence rates of uh ; ph , and uh . k

1

2

3

h

ku  uh k

Order

kp  ph k

Order

ku  uh k

1 4 1 8 1 16 1 32

1.9122e2

–

8.0238e2

–

8.3374e3

4.6805e3

2.0305

1.9913e2

2.0106

1.0202e3

3.0307

1.2166e3

1.9438

5.3908e3

1.8851

1.4301e4

2.8347

3.1021e4

1.9715

1.3984e3

1.9471

1.8834e5

2.9247

1 4 1 8 1 16

3.3183e3

–

1.4504e2

–

1.0009e3

3.8006e4

3.1261

1.7256e3

3.0713

5.5099e5

4.1831

3.8370e5

3.3081

2.2316e4

2.9509

4.1715e6

3.7234

1 4 1 8 1 16

5.1180e4

–

2.8959e3

–

1.3680e4

2.9367e5

4.1233

1.8075e4

4.0019

3.9285e6

5.1219

2.5052e6

3.5512

1.6295e5

3.4715

1.8402e7

4.4160

Order

276

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

0.2

1 0.5

0.1

p

u

0 0

−0.5 −0.1

−1

−0.2 −2

−1.5 −2 −1

−1 0 1 2

y

1

2

−1

0

−2

0 1

x

2

y

1

2

−1

0

−2

x

1

0.18

0.8

0.16

0.6

0.14

0.4

0.12

0.2

0.1

0

0.08

u

p

Fig. 8. Numerical solution of u and p with k ¼ 3; h ¼ 14.

−0.2

0.06

−0.4

0.04

−0.6

0.02

−0.8 −1 0

0

0.5

1 x

1.5

2

−0.02 0

0.5

1 x

1.5

2

0.18

0.18

0.16

0.16

0.14

0.14

0.12

0.12

0.1

0.1

0.08

0.08

u*

u

Fig. 9. Numerical solution of u and p along the line h ¼ p6 with k ¼ 3; h ¼ 14. Also plotted is the exact solution (star).

0.06

0.06

0.04

0.04

0.02

0.02

0

0

−0.02 0

0.5

1 x

1.5

2

−0.02

0

0.5

1 x

1.5

Fig. 10. Numerical solution of u and post processing velocity uh along the line h ¼ p6 with k ¼ 1; h ¼ 14. Also plotted is the exact solution (star).

2

277

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278 Table 6 L2 -errors and convergence rates of uh ; ph ; Lh and uh . k

h

ku  uh k

order

kp  ph k

Order

kL  Lh k

Order

ku  uh k

1.8475e2

–

2.9147e2

–

4.1931e2

–

1.1156e3

1

1 4 1 8 1 16

4.0128e3

2.2029

6.4232e3

2.1820

9.5283e3

2.1377

1.2145e4

3.1994

9.7010e4

2.0484

1.5599e3

2.0418

2.4031e3

1.9873

1.5024e5

3.0150

1 4 1 8 1 16

1.5118e3

–

2.9745e3

–

4.4635e3

–

7.9489e5

1.5617e4

3.2751

3.1959e4

3.2184

4.8770e4

3.1941

4.1198e6

4.2701

1.9219e5

3.0225

4.0598e5

2.9767

6.4532e5

2.9179

2.7019e7

3.9305

1 4 1 8 1 16

1.2126e4

–

3.1428e4

–

4.9263e4

–

5.6349e6

5.7229e6

4.4052

1.5968e5

4.2987

2.7013e5

4.1888

1.5612e7

5.1737

3.7285e7

3.9401

1.0804e6

3.8855

1.9211e6

3.8136

5.5613e9

4.8111

2

1.5

1.5

1

1

0.5

0.5

0

p

u

3

0

−0.5

−0.5

−1

−1

−1.5 −1

−1.5 −1 −0.5

Order

−0.5 0 0.5 1

y

1

0

0.5 x

−0.5

−1

0 0.5 y

1

1

0

0.5

−0.5

−1

x

Fig. 11. Numerical solution of u and p with k ¼ 3; h ¼ 18.

We compare the numerical solution with the exact solution along the line h ¼ p6 in Fig. 9. It shows that the numerical velocity matches very well with the exact one and the jump in pressure is well approximated by the HDG method. The accuracy of the numerical velocity and its postprocessing velocity is compared in Fig. 10. Example 3. In this example, we test a problem in which the interface force has both the normal and tangential components [30]. The exact velocity and pressure are given by

(y r

uðx; yÞ ¼

0 (

v ðx; yÞ ¼ pðx; yÞ ¼

 ry0



r > r0 ;

ð56Þ

r 6 r0 ;

 xr þ rx0

r > r0 ;

0

r 6 r0 ;

ð57Þ

cosðpxÞ cosðpyÞ r > r 0 ; 0

r 6 r0 ;

ð58Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where r ¼ x2 þ y2 . The interface is a circle centered at ð0; 0Þ with the radius r 0 . We choose the computational domain X ¼ ½1; 12 and the radius r0 ¼ 0:5. The external force f and the surface tension can be easily obtained by making the exact solution to satisfy the model problem (1). Obviously, it also have a finite jump across the interface. The L2 -norm of the error of velocity, pressure, gradient of velocity and postprocessing velocity, and convergence rate are presented in Table 6. We can see that optimal convergence rate is obtained for velocity, pressure and the gradient of velocity. The postprocessing velocity has convergence rate of order k þ 2. The numerical solution of u and p with k ¼ 3 and h ¼ 18 is plotted in Fig. 11 (mesh with 1 h ¼ 32 is used in plotting). 6. Concluding remarks In this paper, we investigate the HDG method for solving Stokes interface problems with discontinuous viscosity and variable surface tension. The jump condition on the normal stress ðpI þ mruÞn is naturally incorporated into the HDG

278

B. Wang, B.C. Khoo / Journal of Computational Physics 247 (2013) 262–278

formulation by a constraint on numerical flux accordingly. For polygonal interface, conforming mesh is used. Optimal convergence rates of order k þ 1 are obtained for the velocity, pressure and the gradient of velocity when k-th order polynomials are used to represent the numerical solution. When smooth curved interfaces are involved, we use curvilinear elements to retain the optimal convergence rates. An error estimate for the interface approximation is given based on the CSF model. Numerical results validates that we get a high order and efficient numerical method for solving Stokes interface problem. The investigation of using this method to solve Stokes interface problems with moving interface and multi-phase flows is our on-going work. References [1] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. [2] J. Carrero, B. Cockburn, D. Schötzau, Hybridized, globally divergence-free LDG methods. Part I: The Stokes problem, Math. Comput. 75 (2006) 533–563. [3] Y.C. Chang, T.Y. Hou, B. Merriman, S. Osher, A level set formulation of Eulerian interface capturing methods for incompressible fluid flows, J. Comput. Phys. 124 (1996) 449–464. [4] B. Cockburn, J. Gopalakrishnan, N.C. Nguyen, J. Peraire, F.-J. Sayas, Analysis of HDG methods for Stokes flow, Math. Comput. 80 (2010) 723–760. [5] B. Cockburn, J. Gopalakrishnan, The derivation of hybridizable discontinuous Galerkin methods for Stokes flow, SIAM J. Numer. Anal. 47 (2009) 1092– 1125. [6] B. Cockburn, J. Gopalakrishnan, R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009) 1319–1365. [7] B. Cockburn, G. Kanschat, D. Schötzau, C. Schwab, Local discontinuous Galerkin methods for the Stokes system, SIAM J. Numer. Anal. 40 (2002) 319– 343. [8] R. Cortez, The method of regularized Stokeslets, SIAM J. Sci. Comput. 23 (2001) 1204–1225. [9] S. Gross, A. Reusken, Finite elment discretization error analysis of a surface tension force in two-phase incompressible flows, SIAM J. Numer. Anal. 45 (2007) 1679–1700. [10] S. Gross, A. Reusken, Numerical Methods for Two-phase Incompressible Flows, Springer Press, New York, USA, 2011. [11] J.S. Hesthaven, T. Warburton, Nodal Discontinuous Galerkin Methods Algorithms, Analysis, and Applications, Springer, 2008. ISBN: 978-0-387-720644. [12] G. Hou, J. Wang, A. Layton, Numerical methods for fluid–structure interaction— a review, Commun. Comput. Phys. 12 (2012) 337–377. [13] L.N.T. Huynh, N.C. Nguyen, J. Paraire, B.C. Khoo, A high order hybridizable discontinuous Galerkin method for elliptic interface problems, Int. J. Number. Meth. Eng. 93 (2013) 183–200. [14] M-C. Lai, H-C. Tseng, A simple implementation of the immersed interface methods for Stokes flows with singular forces, Comput. Fluids 37 (2008) 99– 106. [15] D.V. Le, B.C. Khoo, J. Peraire, An immersed interface method for viscous incompressible flows involving rigid and flexible boundaries, J. Comput. Phys. 220 (2006) 109–138. [16] R.J. LeVeque, Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal. 31 (1994) 1019–1044. [17] R.J. LeVeque, Z. Li, Immersed interface method for Stokes flow with elastic boundaries or surface tension, SIAM J. Sci. Comput. 18 (1997) 709–735. [18] L. Li, RJ. LeVeque, An immersed interface method for incompressible Navier–Stokes equations, SIAM J. Sci. Comput. 25 (2003) 832–856. [19] Z. Li, M.C. Lai, The immersed interafce method for the Navier–Stokes equations with singular forces, J. Comput. Phys. 171 (2001) 822–842. [20] S. Nagrath, K.E. Jansen, R.T. Lahey, Computation of incompressible bubble dynamics with a stabilized finite element level set method, Comput. Meth. Appl. Mech. Eng. 194 (2005) 4565–4587. [21] N.C. Nguyen, J. Peraire, B. CockBurn, Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell’s equations, J. Comput. Phys. 230 (2011) 7151–7175. [22] N.C. Nguyen, J. Peraire, B. CockBurn, An implicit high-order hybridizable discontinuous Galerkin method for the incompressible Navier–Stokes equations, J. Comput. Phys. 230 (2011) 1147–1170. [23] N.C. Nguyen, J. Peraire, B. CockBurn, A hybridizable discontinuous Galerkin method for Stokes flow, Comput. Meth. Appl. Mech. Eng. 199 (2010) 582– 597. [24] N.C. Nguyen, J. Peraire, B. CockBurn, A hybridizable discontinuous Galerkin method for the incompressible Navier–Stokes equations (AIAA Paper 2010362), in: Proceedings of the 48th AIAA Aerospace Sciences Meeting and Exhibit, Orlando, FL, Janury 2010. [25] E. Olsson, G. Kreiss, A conservative level set method for two phase flow, J. Comput. Phys. 210 (2005) 225–246. [26] M.A. Olshanskii, A. Reusken, Analysis of a Stokes interface problem, Numer. Math. 103 (2006) 129–149. [27] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (1977) 220–252. [28] C.S. Peskin, The immersed boundary method, Acta Numer. 11 (2002) 479–517. [29] A.M. Roma, C.S. Peskin, M.J. Berger, An adaptive version of the immersed boundary method, J. Comput. Phys. 153 (1999) 509–534. [30] Z.J. Tan, K.M. Lim, B.C. Khoo, A fast immersed interface method for solving Stokes flows on irregular domains, Comput. Fluids 38 (2009) 1973–1983. [31] Z.J. Tan, K.M. Lim, B.C. Khoo, An immersed interface method for Stokes flows with fixed/moving interfaces and rigid boundaries, J. Comput. Phys. 228 (2009) 6855–6881.