A staggered semi-implicit discontinuous Galerkin method for the two dimensional incompressible Navier–Stokes equations

Applied Mathematics and Computation 248 (2014) 70–92

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A staggered semi-implicit discontinuous Galerkin method for the two dimensional incompressible Navier–Stokes equations Maurizio Tavelli a, Michael Dumbser b,⇑ a b

Department of Mathematics, University of Trento, Via Sommarive 14, I-38050 Trento, Italy Laboratory of Applied Mathematics, Department of Civil, Environmental and Mechanical Engineering, University of Trento, Via Mesiano 77, I-38123 Trento, Italy

a r t i c l e

i n f o

Keywords: Semi-implicit discontinuous Galerkin schemes Staggered unstructured triangular meshes High order staggered finite element schemes Curved isoparametric elements Incompressible Navier–Stokes equations

a b s t r a c t In this paper we propose a new spatially high order accurate semi-implicit discontinuous Galerkin (DG) method for the solution of the two dimensional incompressible Navier–Stokes equations on staggered unstructured curved meshes. While the discrete pressure is defined on the primal grid, the discrete velocity vector field is defined on an edge-based dual grid. The flexibility of high order DG methods on curved unstructured meshes allows to discretize even complex physical domains on rather coarse grids. Formal substitution of the discrete momentum equation into the discrete continuity equation yields one sparse linear equation system with four non-zero blocks per element for only one scalar unknown, namely the pressure. The method is computationally efficient, since the resulting system is not only very sparse but also symmetric and positive definite for appropriate boundary conditions. Furthermore, all the volume and surface integrals needed by the scheme presented in this paper depend only on the geometry and the polynomial degree of the basis and test functions and can therefore be precomputed and stored in a preprocessor stage, which leads to savings in terms of computational effort for the time evolution part. In this way also the extension to a fully curved isoparametric approach becomes natural and affects only the preprocessing step. The method is validated for polynomial degrees up to p ¼ 3 by solving some typical numerical test problems and comparing the numerical results with available analytical solutions or other numerical and experimental reference data. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction The main difficulty in the numerical solution of the incompressible Navier–Stokes equations lies in the pressure Poisson equation and the associated linear equation system to be solved on the discrete level. This is closely related to the elliptic nature of these equations, where boundary conditions affect instantly the solution everywhere inside the domain. While finite difference schemes for the incompressible Navier–Stokes equations are well-established for several decades now [50,63,62,78], as well as continuous finite element methods [73,11,55,44,79,52,53], the development of high order discontinuous Galerkin (DG) finite element methods for the incompressible Navier–Stokes equations is still a very active topic of ongoing research.

⇑ Corresponding author. E-mail addresses: [email protected] (M. Tavelli), [email protected] (M. Dumbser). http://dx.doi.org/10.1016/j.amc.2014.09.089 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

71

Several high order DG methods for the incompressible Navier–Stokes equations have been recently presented in literature, see for example [4,69,43,61,65,66,34,58], or the work of Bassi et al. [3] based on the technique of artificial compressibility, originally introduced by Chorin [24,25]. In this paper we propose a new, spatially high order accurate semi-implicit DG finite element scheme that is based on the general ideas of [39,71], following the philosophy of semi-implicit staggered finite difference schemes, which have been successfully used in the past for the solution of the incompressible Navier–Stokes equations [50,63,62,78] and the free surface shallow water and Navier–Stokes equations, see [54,18,19,21,80,15]. Very recent developments in the field of such semi-implicit finite difference schemes for the free surface Navier–Stokes equations can be found in [16,20,17], together with their theoretical analysis presented in [12–14,23]. In our semi-implicit staggered DG scheme, the discrete pressure is defined on the control volumes of the primal triangular mesh, while the discrete velocity vector is defined on an edge-based, quadrilateral dual mesh. The nonlinear convective terms are discretized explicitly in time, using a classical RKDG scheme [32,31,33] based on the local Lax-Friedrichs (Rusanov) flux [67], while the viscous terms are discretized implicitly using a fractional step method. The DG discretization of the viscous fluxes is based on the formulation of Gassner et al. [46], who obtained the viscous numerical flux from the solution of the Generalized Riemann Problem (GRP) of the diffusion equation. The solution of the GRP has first been used to construct numerical methods for hyperbolic conservation laws by Ben Artzi and Falcovich [6] and Toro and Titarev [76,74]. The discrete momentum equation is then inserted into the discrete continuity equation in order to obtain the discrete form of the pressure Poisson equation. The chosen dual grid used here is taken as the one used in [7,75,26,71,8], which leads to a sparse linear system for the scalar pressure with only four non-zero blocks per element. Once the new pressure field is known, the velocity vector field can subsequently be updated directly. Very recently, an accurate and efficient pressurebased hybrid finite volume/ finite element solver using staggered unstructured meshes has been proposed in [8]. Other staggered DG schemes have been used in [26,27,30,28,29,59,60]. However, to our knowledge, none of these schemes has ever been applied to the incompressible Navier–Stokes equations. To our knowledge, a staggered DG scheme has been proposed only for the Stokes system so far, see [57]. For alternative semi-implicit DG schemes on collocated grids see [77,47,35–37]. The rest of the paper is organized as follows: in Section 2 the numerical method is described in detail, while in Section 3 a set of numerical test problems is solved in order to study the accuracy of the presented approach. Some concluding remarks are given in Section 4. 2. DG scheme for the 2D incompressible Navier–Stokes equations 2.1. Governing equations The two dimensional incompressible Navier–Stokes equations and the continuity equation are given by

@v þ r  Fc þ rp ¼ mDv ; @t r  v ¼ 0;

ð1Þ ð2Þ

where p ¼ P=q indicates the normalized fluid pressure; P is the physical pressure and q is the constant fluid density; m is the kinematic viscosity coefficient; v ¼ ðu; v Þ is the velocity vector, where u and v are the velocity components in the x and y direction, respectively; Fc ¼ v  v is the flux tensor of the nonlinear convective terms, namely:

 Fc ¼

uu

uv

v u vv

 :

The viscosity term is first written as mDv ¼ r  ðmrv Þ and then grouped with the nonlinear convective term. So Eq. (1) becomes

@v þ r  F þ rp ¼ 0; @t

ð3Þ

where F ¼ Fðv ; rv Þ ¼ Fc ðv Þ  mrv is a nonlinear tensor that depends on the velocity and its gradient, see e.g. [46,38]. We further use the abbreviation Lðv Þ ¼ @t@ v þ r  F. 2.2. Unstructured grid In this paper we use the same general unstructured staggered mesh proposed in [71]. In this section we briefly summarize the grid construction and the main notation. The computational domain is covered with a set of N i non-overlapping triangles T i with i ¼ 1 . . . N i . By denoting with N j the total number of edges, the jth edge will be called Cj . BðXÞ denotes the set of indices j corresponding to boundary edges. The three edges of each triangle T i constitute the set Si defined by Si ¼ fj 2 ½1; N j  j Cj is an edge of T i g. For every j 2 ½1 . . . N j   BðXÞ there exist two triangles i1 and i2 that share Cj . It is possible to assign arbitrarily a left and a right triangle called ‘ðjÞ and rðjÞ, respectively. The standard positive direction is

72

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

assumed to be from left to right. Let ~ nj denote the unit normal vector defined on the edge j and oriented with respect to the positive direction from left to right. For every triangular element i and edge j 2 Si , the neighbor triangle of element T i at edge Cj is denoted by }ði; jÞ. For every j 2 ½1; Nj   BðXÞ the quadrilateral element associated to j is called Rj and it is defined, in general, by the two centers of gravity of ‘ðjÞ and rðjÞ and the two terminal nodes of Cj , see also [7,75,71]. We denote by T i;j ¼ Rj \ T i the intersection element for every i and j 2 Si . Fig. 1 summarizes the notation used here, the main triangular and the dual quadrilateral meshes. According to [71], we call the mesh of triangular elements fT i gi2½1;Ni  the main grid or the primal grid and the quadrilateral grid fRj gj2½1;Nj  is called the dual grid. The dual grid is obtained by joining the vertices of each edge of the primal mesh with two points inside the left and right triangles adjacent to this edge, see Fig. 1. Typically the internal point is chosen as the triangle barycenter, but also the center of the circumcircle of each primary triangle could be chosen, which leads to an orthogonal dual grid. On the dual grid we define the same quantities as for the main grid, briefly: N l is the total amount of edges of Rj ; Cl indicates the lth edge; 8j, the set of edges l of j is indicated with Sj ; 8l; ‘jl ðlÞ and rjl ðlÞ are the left and the right quadrilateral element, respectively; 8l; ~ nl is the standard normal vector defined on l and assumed positive with respect to the standard orientation on l (defined, as above, from the left to the right). Finally, each triangle T i is defined starting from an arbitrary node and oriented in counter-clockwise direction. Similarly, each quadrilateral element Rj is defined, starting from ‘ðjÞ and oriented in counter-clockwise direction. 2.3. Basis functions According to [71] we proceed as follows: we first construct the polynomial basis up to a generic polynomial degree p on some reference triangular and quadrilateral elements. In order to do this we take T std ¼ fðn; cÞ 2 R2;þ j c 6 1  n _ 0 6 n 6 1g as the reference triangle and the unit square as the reference quadrilateral element Rstd ¼ ½0; 12 . Using the standard nodal approach of conforming continuous finite elements, we obtain N / ¼ ðpþ1Þðpþ2Þ basis functions f/k gk2½1;N/  on T std and 2 N w ¼ ðp þ 1Þ2 basis functions on Rstd . The connection between reference and physical space is performed by the maps T 1 i

Tj

Ti

T 1 j

T i ! T std for every i ¼ 1 . . . N i ; Rj ! Rstd for every j ¼ 1 . . . N j and its inverse, called T i  T std and Rj  Rstd , respectively. The maps from physical coordinates to reference coordinates can be constructed following a classical sub-parametric or a complete iso-parametric approach. 2.4. Semi-implicit DG scheme We define the spaces of piecewise polynomials used on the main grid and the dual grid as follows, p Vm h ¼ f/ : /jT i 2 P ðT i Þ; 8i 2 ½1; N i g;

and V dh ¼ fw : wjRj 2 Qp ðRj Þ;

8j 2 ½1; Nj   BðXÞg;

ð4Þ

Fig. 1. Example of a triangular mesh element with its three neighbors and the associated staggered edge-based dual control volumes, together with the notation used throughout the paper.

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

73

where Pp ðT i Þ is the space of polynomials of degree at most p on T i , while Qp ðRj Þ is the space of tensor products of one-dimensional polynomials of degree at most p on Rj . The discrete pressure ph is defined on the main grid while the discrete velocity vector field v h is defined on the dual grid, d namely ph 2 V m h and v h 2 V h for each component of the velocity vector. The numerical solution of (2)–(3) is represented by piecewise polynomials and written in terms of the basis functions on the primary and the dual grid as

pi ðx; y; tÞ ¼

N/ X ðiÞ ^ i ðtÞ; ^l;i ðtÞ ¼: /ðiÞ ðx; yÞp /l ðx; yÞp

ð5Þ

l¼1

vj ðx; y; tÞ ¼

Nw X ðjÞ ^ l;j ðtÞ ¼: wðjÞ ðx; yÞv ^ j ðtÞ; wl ðx; yÞv

ð6Þ

l¼1

where the vector of basis functions /ðx; yÞ and wðx; yÞ are generated from /ðn; cÞ on T std , and from wðn; cÞ on Rstd , respectively. Formally /ðiÞ ðx; yÞ ¼ /ðT i ðx; yÞÞ for i ¼ 1 . . . N i and wðjÞ ðx; yÞ ¼ wðT j ðx; yÞÞ for every j ¼ 1 . . . N j . A weak formulation of Eq. (2) is obtained by multiplying it by / and integrating over a control volume T i , for every k ¼ 1 . . . N / . The resulting weak formulation of (2) reads

Z

Ti

ðiÞ

/k r  v dx ¼ 0

ð7Þ

with dx ¼ dxdy. Similarly, multiplication of the momentum Eq. (3) by w and integrating over a control volume Rj one obtains, componentwise,

Z Rj



ðjÞ

wk

 Z @v ðjÞ þ r  F dx þ wk rp dx ¼ 0 @t Rj

ð8Þ

for every j ¼ 1 . . . N j and k ¼ 1 . . . N w . Using integration by parts Eq. (7) yields

I

ðiÞ

@T i

/k

v  ~ni ds 

Z

Ti

r/ðiÞ k  v dx ¼ 0;

ð9Þ

ni indicates the outward pointing unit normal vector. The discrete pressure ph in general presents a discontinuity on where ~ v h jumps on the edges of Rj (see Fig. 2). Hence, Eqs. (8) and (9) have to be split as follows:

Cj and also the discrete velocity field X Z

ðiÞ

Cj

j2Si

/k

v j  ~nij ds 

!

Z T i;j

r/kðiÞ  v j dx ¼ 0

ð10Þ

and

Z Rj



ðjÞ

wk

 Z Z Z   @ vj ðjÞ ðjÞ ðjÞ þ r  Fj dx þ wk rp‘ðjÞ dx þ wk rprðjÞ dx þ wk prðjÞ  p‘ðjÞ ~ nj ds ¼ 0; @t T ‘ðjÞ;j T rðjÞ;j Cj

where ~ nij ¼ ~ ni jCj . Definitions (5) and (6) allow to rewrite the above equations by splitting the spatial and temporal variables, namely

X Z j2Si

Cj

ðiÞ ðjÞ /k wl ~ nij ds

v^ l;j 

Z T i;j

!

r

ðiÞ ðjÞ /k wl dx

v^ l;j

¼0

Fig. 2. Jumps of p on the main grid (left) and of ~ v on the dual grid (right).

ð11Þ

74

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

and

Z Rj

ðjÞ ðjÞ ^ l;j Þ þ wk wl dxLh ðv

Z

ðjÞ

T ‘ðjÞ;j

Z

^l;‘ðjÞ þ dx p

ðjÞ

~ ^l;rðjÞ  nj ds p

Z

þ

ð‘ðjÞÞ

wk r/l

Cj

ðjÞ

T rðjÞ;j

ðrðjÞÞ

wk /l

Z

Cj

ðrðjÞÞ

wk r/l ðjÞ

^l;rðjÞ dx p

ð‘ðjÞÞ

wk /l

~ ^l;‘ðjÞ ¼ 0; nj ds p

ð12Þ

where we have used the standard summation convention for the repeated index l. Lh is an appropriate discretization of the operator L and will be given later. For every i and j, Eqs. (11)–(12) are written in a compact matrix form such as

X ^j ¼ 0 Di;j v

ð13Þ

j2Si

and

^ rðjÞ  Lj p ^ ‘ðjÞ ¼ 0; ^ j Þ þ Rj p M j Lh ðv

ð14Þ

respectively, where:

Mj ¼

Z

ðjÞ

Rj

Z

Di;j ¼

ðiÞ ðjÞ /k wl ~ nij ds 

Cj

Rj ¼

Z

ðjÞ

Cj

Lj ¼

ðjÞ

wk wl dx;

ðjÞ

Cj

Z

ðrðjÞÞ

~ nj ds þ

wk /l

Z

ð15Þ

T i;j

r/kðiÞ wðjÞ l dx;

Z

ðjÞ

T rðjÞ;j

ð‘ðjÞÞ

~ nj ds 

wk /l

Z

ðrðjÞÞ

wk r/l ðjÞ

T ‘ðjÞ;j

ð16Þ

ð‘ðjÞÞ

wk r/l

dx;

ð17Þ

dx:

ð18Þ

According to [71] the action of tensors L and R can be generalized by introducing the new tensor Qi;j , defined as

Qi;j ¼

Z

T i;j

where

ðjÞ

ðiÞ

wk r/l dx 

Z

Cj

ðjÞ

ðiÞ

wk /l

ri;j~ nj ds;

ð19Þ

ri;j is a sign function defined by

ri;j ¼

rðjÞ  2i þ ‘ðjÞ : rðjÞ  ‘ðjÞ

ð20Þ

In this way Q‘ðjÞ;j ¼ Lj and QrðjÞ;j ¼ Rj , and then Eq. (14) becomes in terms of Q

^ rðjÞ þ Q‘ðjÞ;j p ^ ‘ðjÞ ¼ 0; ^ j Þ þ QrðjÞ;j p M j Lh ðv

ð21Þ

or, equivalently,

^ i þ Q}ði;jÞ;j p ^ }ði;jÞ ¼ 0: ^ j Þ þ Qi;j p M j Lh ðv

ð22Þ

We discretize the velocity in Eq. (13) implicitly and the pressure in Eq. (14) semi-implicitly by using the theta method in time, namely

8P ^ nþ1 ¼ 0; < j2Si Di;j v j nþ1 b n ^ v  F v : ^ nþh ^ nþh M j j Dt j þ QrðjÞ;j p rðjÞ þ Q‘ðjÞ;j p‘ðjÞ ¼ 0; ^ nþh ¼ hp ^ nþ1 þ ð1  hÞp ^ n ; and h is an implicitness factor to be taken in the range h 2 where p Eqs. (23) as described above and using the formulation (22), we get for every i and j 2 Si

X ^ nþ1 Di;j v ¼ 0; j

ð23Þ 1 2

 ; 1 , see e.g. [18]. Discretizing

ð24Þ

j2Si

Mj

v^ nþ1  Fc v nj j Dt

    ^ nþ1 ^ nþ1 ^n ^n þ h Qi;j p þ Q}ði;jÞ;j p i }ði;jÞ þ ð1  hÞ Qi;j pi þ Q}ði;jÞ;j p}ði;jÞ ¼ 0;

ð25Þ

where Fc v nj is an appropriate discretization of the nonlinear convective and viscous terms. The details for the computation of n c F v j will be presented later in Section 2.5. Formal substitution of the momentum equation (25) into the continuity equation (24), see also [19,39], yields

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

hDt

X X n ^ nþ1  hDt Di;j M 1 ^ nþ1 Di;j M 1 j Qi;j pi j Q}ði;jÞ;j p}ði;jÞ ¼ bi ; j2Si

75

ð26Þ

j2Si

where

 X X   n bi ¼  Di;j Fc v nj þ ð1  hÞDt Di;j M j 1 Qi;j p^ ni þ Q}ði;jÞ;j p^ n}ði;jÞ;j ; j2Si

ð27Þ

j2Si

groups all the known terms at time tn . ^ nþ1 Eq. (26) represents a block four-point system for the new pressure p . It can be interpreted as the discrete form of the i pressure Poisson equation of the incompressible Navier–Stokes equations. Once the new pressure field is known, the velocity field can be readily updated from the momentum equation, Eq. (25). We emphasize that in the present algorithm, the only unknown is the scalar pressure ph . It remains to complete the system by introducing the boundary conditions. In order to do this observe how, for i 2 ½1; N i  and j 2 Si \ BðXÞ, the boundary element Rj ¼ T i;j is a triangular element and not a quadrilateral element. The basis functions to be used are the one generated on T std . In this way the matrices M j ; Di;j ; Qi;j defined in (15), (16) and (19), have to be modified for boundary elements. For every j 2 Si \ BðXÞ

vj ¼

N/ X ^ l;j ; /l v

ð28Þ

l¼1

where the /l are the basis functions on the reference triangle T std . The matrices can be recomputed for j 2 Si \ BðXÞ and will be called D@i;j ; Q@i;j . Eqs. (26) and (27) are consequently computed with the triangular boundary elements and one so obtains

"

X

hDt 

j2Si \BðXÞ

#

X

@ D@i;j M 1 j Qi;j 

^ nþ1  hDt Di;j M 1 j Qi;j pi

j2Si BðXÞ

X

~n ^ nþ1 Di;j M 1 j Q}ði;jÞ;j p}ði;jÞ ¼ bi ;

ð29Þ

j2Si BðXÞ

where now the vector of known terms is

X

~n ¼  b i

Di;j Fc v nj þ

j2Si BðXÞ

D@i;j Fc v nj

j2Si \BðXÞ

X

þ ð1  hÞDt

X

X

@ ^n D@i;j M 1 j Qi;j pi þ ð1  hÞDt

j2Si \BðXÞ

  ^ ni þ Q}ði;jÞ;j p ^ n}ði;jÞ : Di;j M 1 Qi;j p j

ð30Þ

j2Si BðXÞ

As implied by Eq. (29), the stencil of the present scheme only involves the ith element and its direct Neumann neighbors. Thus, since #Si ¼ 3, the system described by (29) is a sparse system with four non-zero blocks per element. As we will show later, the system is symmetric and positive definite for appropriate boundary conditions, hence it can be efficiently solved by using a matrix-free implementation of the conjugate gradient algorithm [51]. Once the new pressure has been computed, the new velocity field can be readily updated from Eq. (25) for every j R BðXÞ:

v^ nþ1 ¼ Fc v nj  hDtM 1 j j

    1 ^ nþ1 ^ nþ1 ^ ni þ Q}ði;jÞ;j p ^ n}ði;jÞ : Qi;j p þ Q}ði;jÞ;j p Qi;j p i }ði;jÞ  ð1  hÞDtM j

ð31Þ

The above Eqs. (29)–(31) can be modified for j 2 BðXÞ according to the type of boundary conditions (velocity or pressure boundary condition). Note that all the matrices used in the above algorithm can be precomputed once and forall for a given mesh and polynomial degree p. 2.5. Nonlinear convection–diffusion In problems where the convective term and the viscosity can be neglected we can take Fc v nj ¼ v^ nj in Eq. (27). Otherwise, an explicit cell-centered RKDG method [32] on the dual mesh is used in this paper for the discretization of the nonlinear convective terms. The viscosity contribution is discretized implicitly using a fractional step method, in order to avoid additional restrictions on the time step Dt. The semi-discrete DG scheme for the nonlinear convection–diffusion terms on the dual mesh is given by

Z

wk Rj

d v h dx þ dt

Z

n ds  wk Gh  ~

@Rj

Z

rwk  Fðv h ; rv h Þ dx ¼ 0

ð32Þ

Rj

and the numerical flux for both, the convective and the viscous contribution, is given by [67,46,38] as

Gh  ~ n¼

 1 1 Fðv hþ ; rv hþ Þ þ Fðv h ; rv h Þ  ~ n  smax v hþ  v h 2 2

ð33Þ

with

 smax ¼ 2 max jv h  ~ nj; jv hþ  ~ nj þ

2m 2p þ 1 pffiffipffi ; þ  h þh 2

ð34Þ

76

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

which contains the maximum eigenvalue of the Jacobian matrix of the purely convective transport operator Fc in normal  direction, see [39], and the stabilization term for the viscous flux, see [38,46]. Furthermore, the v  h and rv h denote the velocity vectors and their gradients, extrapolated to the boundary of Rj from within the element Rj and from the neighbor þ  element, respectively. h and h are the maximum radii of the inscribed circle in Rj and the neighbor element, respectively. A classical third order accurate TVD Runge–Kutta method is used for time integration of the nonlinear convective terms, see e.g. [70,49,32], since the explicit discretization of higher order DG schemes with a simple first order Euler method in time would lead to a linearly unstable scheme. The above method requires that the time step size is restricted by a CFL-type restriction for DG schemes, namely:

Dt ¼

CFL hmin  ; 2p þ 1 2jv max j

ð35Þ

where hmin is the smallest incircle diameter; CFL < 0.5; and v max is the maximum convective speed. Furthermore, the time step of the global semi-implicit scheme is not affected by the local time step used for the time integration of the convective terms if a local time stepping/subcycling approach is employed, see [22,72]. Implicit discretization of the viscous terms in (32) with a fractional step method involves a block five-point system that can be efficiently solved using the GMRES algorithm [68]. The solution of this system is not necessary in problems where the viscous terms are small and can be integrated explicitly in time. The stability of the method is linked to the nonlinear convective term, so the method is stable under condition (35). 2.6. Extension to curved elements The maps used to switch between reference and physical space can be defined using a simple sub-parametric vertex based approach or a fully isoparametric approach. In the first case only the vertices of the elements T i and Rj are required to map the physical element into the reference one and vice versa. In this simple case an explicit expression for the maps T i ; T 1 and T 1 can be computed while for the map T j we use the Newton method (see e.g. [71]). A simple extension to i j the complete isoparametric case requires to store more information about each element, namely we need to know the coordinates of the nodes fðX; YÞik gk¼1;N/ for each triangular element i and fðX; YÞjl gl¼1;Nw for each quadrilateral element j. The inverse maps T 1 and T 1 are defined by using the same basis functions /k and wk used for representing the discrete i j solution of the PDE, i.e. we have

x¼

N/ X

/k X ik ;

y¼

k

N/ X

/k Y ik

ð36Þ

wk Y jk

ð37Þ

k

and

x¼

Nw X

wk X jk ;

y¼

k

Nw X k

for triangles and quadrilateral elements, respectively. In this case the maps T i and T j become nonlinear and so the Newton method has to be used for both. Also the Jacobian and the normal vectors are not, in general, constant through the element and the edges, respectively. The main advantage of this approach is that now the edges become curved and so the computational domain can better approximate the physical one. It is important to observe how this approach affects only the preprocessing step. 2.7. Remarks on the main system and further improvements In this section we will show how the main system for the computation of the pressure, developed in Section 2.4 results symmetric and, in general, positive semi-definite. These results allow us to use very fast iterative methods to solve the linear system, such as the conjugate gradient method [51], with a significant gain in terms of computational time. In order to do this observe how, from the definitions (16) and (19), we can further generalize the action of D in terms of Q such as D ¼ Q> since

Q>i;j

Z

¼ ¼

Z

Xi;j

Xi;j

ðjÞ wk ðjÞ

r

ðiÞ /l dx ðiÞ



wl r/k dx þ

!>

Z

Z

Cj

Cj

ðjÞ ðiÞ wk /l ðjÞ

ðiÞ

wl /k

ri;j~ nj ds

ri;j~ nj ds ¼ Di;j

nj , else, it is ~ nj ; 8i; j 2 Si . Consequently, the main system (26) can be written as and if i ¼ ‘ðjÞ; ~ nij coincides with ~

A : hDt

X >  1 X  1 n ^ nþ1 ^ nþ1 Qi;j M j Qi;j p þ hDt Q>i;j M j Q}ði;jÞ;j p i }ði;jÞ ¼ bi ; j2Si

j2Si

ð38Þ

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

77

that we will call in the following A. If we do not introduce any boundary conditions, we have the following Lemma 1. Without any boundary conditions the system A is singular. Proof 1. Let ph 2 V m h , in order to show that A is singular we investigate the kernel of the linear operator A. Since detA – 0 () KerA ¼ f0g, we would like to show that the kernel does not contain only the zero. A weak formulation of rph over Xj is given by Q‘ðjÞ;j p‘ðjÞ þ QrðjÞ;j prðjÞ , then we have the identity

^ ‘ðjÞ þ QrðjÞ;j p ^ rðjÞ  0 () rpjX ¼ 0: Q‘ðjÞ;j p j

ð39Þ

We are looking for a ph – 0 such that Aph ¼ 0. For a fixed i 2 ½1; N i ,

 h Dt

X j2Si

 h Dt

X j2Si

 h Dt

X

X  1  1 ^ nþ1 ^ nþ1 Di;j M j Qi;j p  hDt Di;j M j Q}ði;jÞ;j p i }ði;jÞ ¼ 0; j2Si

i  1 h ^ nþ1 ^ nþ1 Di;j M j Qi;j p þ Q}ði;jÞ;j p i }ði;jÞ ¼ 0; i  1 h ^ nþ1 ^ nþ1 ¼ 0: Di;j M j Q‘ðjÞ;j p ‘ðjÞ þ QrðjÞ;j prðjÞ

ð40Þ

j2Si

Hence, if p ¼ constant, the left side of (40) vanishes and then fpi  c 8i; c 2 Rg  kerA. h This represents a natural result since the incompressible Navier-Stokes equations depend only on the gradient of the pressure and not directly on the pressure. Once we have an exact solution for the pressure pe , then every solution of the kind pe þ c with c 2 R is also a solution. If we introduce the boundary conditions and we specify the pressure in at least one point (i.e. in at least one degree of freedom), this is equivalent to choose the constant c and the system becomes non-singular. The following results state that the developed system has several important properties such as the symmetry and, in general, positive semi-definiteness. Lemma 2 (Symmetry). The system matrix of A is symmetric. Proof 2. In the following we denote with ði; kÞ the kth degree of freedom of the ith element. For the symmetry of A we have ~ as ð~i; kÞ ~ act on ði; kÞ. If i ¼ ~i, the action is described by P Q> M j 1 Q that is trivially to verify that ði; kÞ act on ð~i; kÞ i;j i;j j2Si ~ symmetric since M j ¼ M > j is symmetric. If i R }ði; Si Þ the two actions are zero so it is also trivially verified. Remains the 1 case ~i 2 }ði; Si Þ. In this case, the actions of the right element on the left one and vice versa are, respectively, Q> ‘ðjÞ;j M j QrðjÞ;j 1 > and QrðjÞ;j M j Q‘ðjÞ;j . A simple computation leads to

M 1 j QrðjÞ;j ðk; lÞ ¼

Nw X M 1 8k ¼ 1 . . . N w ; l ¼ 1 . . . N / j ðk; nÞQrðjÞ;j ðn; lÞ n¼1

and then 8k ¼ 1 . . . N / ; l ¼ 1 . . . N / ,

Q>‘ðjÞ;j M 1 j QrðjÞ;j ðk; lÞ ¼

Nw   X Q>‘ðjÞ;j ðk; cÞ M 1 j QrðjÞ;j ðc; lÞ

c¼1

¼

Nw Nw X X Q>‘ðjÞ;j ðk; cÞ M 1 j ðc; nÞQrðjÞ;j ðn; lÞ

c¼1

¼

Nw X

n¼1

Q>‘ðjÞ;j ðk; cÞM 1 j ðc; nÞQrðjÞ;j ðn; lÞ

c;n¼1

¼

Nw X

> Q‘ðjÞ;j ðc; kÞM 1 j ðc; nÞQrðjÞ;j ðl; nÞ

c;n¼1

¼

Nw X

Q>rðjÞ;j ðl; nÞM 1 j ðn; cÞQ‘ðjÞ;j ðc; kÞ

c;n¼1

¼ Q>rðjÞ;j M 1 j Q‘ðjÞ;j ðk; lÞ:



Lemma 3. The matrix A is positive semi-definite, i.e. x> Ax P 0 8x 2 RNi N/ .

ð41Þ

78

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

Proof 3. We do the computation directly. x> Ax ¼

P

i ðx

>

AxÞi and

X X ðx> AxÞi ¼ xi Q>i;j M 1 Q>i;j M 1 j Qi;j xi þ xi j Q}ði;jÞ;j x}ði;jÞ j2Si

j2Si

>  1  X 1 ¼ M j 2 Qi;j xi M j 2 Qi;j xi j2Si

>  1  X 1 M j 2 Qi;j xi M j 2 Q}ði;jÞ;j x}ði;jÞ ;

þ

j2Si

where we used that M j is symmetric and positive definite, hence M 1 is symmetric and positive definite and then exists the j  1 >  1  12 1 1 so called matrix square root operator, namely 9M j such that M j ¼ M j 2 M j 2 . By defining T i;j :¼ M j 2 Qi;j we obtain

ðx> AxÞi ¼

X

T i;j xi

> 

X >  T i;j xi þ T i;j xi T }ði;jÞ;j x}ði;jÞ :

j2Si

ð42Þ

j2Si

and consequently Ni X X 

x> Ax ¼

T i;j xi

Ni X >  X  >  T i;j xi þ T i;j xi T }ði;jÞ;j x}ði;jÞ :

i¼1 j2Si

ð43Þ

i¼1 j2Si

Remark that the double summation

PN i P i¼1

j2Si

sum every element i and edge j. From the edge point of view, every edge

gives two contributions, one given when i ¼ ‘ðjÞ and one when i ¼ rðjÞ. The double summation can be consequently inverted as follows: Ni X X 

T i;j xi

> 

T i;j xi

> 

N

N

j j X  >  X  >  T i;j xi ¼ T ‘ðjÞ;j x‘ðjÞ T ‘ðjÞ;j x‘ðjÞ þ T rðjÞ;j xrðjÞ T rðjÞ;j xrðjÞ

i¼1 j2Si

j¼1

Ni X X 

T }ði;jÞ;j x}ði;jÞ



j¼1 Nj Nj X  >  X  >  T ‘ðjÞ;j x‘ðjÞ T rðjÞ;j xrðjÞ þ T rðjÞ;j xrðjÞ T ‘ðjÞ;j x‘ðjÞ ¼

i¼1 j2Si

j¼1

ð44Þ

j¼1

and then, by recomposing everything

x> Ax ¼

Nj h X  >   >   >   >  i T ‘ðjÞ;j x‘ðjÞ T ‘ðjÞ;j x‘ðjÞ þ T rðjÞ;j xrðjÞ T rðjÞ;j xrðjÞ þ T ‘ðjÞ;j x‘ðjÞ T rðjÞ;j xrðjÞ þ T rðjÞ;j xrðjÞ T ‘ðjÞ;j x‘ðjÞ j¼1

¼

Nj X 

T ‘ðjÞ;j x‘ðjÞ þ T rðjÞ;j xrðjÞ

>  T ‘ðjÞ;j x‘ðjÞ þ T rðjÞ;j xrðjÞ

j¼1

¼

" Nj T ‘ðjÞ;j X 0

j¼1

¼

¼

!

0



T rðjÞ;j

Nj T ‘ðjÞ;j X ðx‘ðjÞ ; xrðjÞ Þ 0 j¼1

x‘ðjÞ

!#> "

T ‘ðjÞ;j

0

!

x‘ðjÞ

!#

 xrðjÞ xrðjÞ 0 T rðjÞ;j !> ! ! x‘ðjÞ 0 T ‘ðjÞ;j 0 T rðjÞ;j

0

T rðjÞ;j

xrðjÞ

Nj X ~ x>j T > T ~ xj :

ð45Þ

j¼1 > ~ :¼ T > T is a positive semi-definite matrix by construction, ~ ~~ And, since T x> j T xj P 0 and then x Ax ¼

We introduce now the boundary elements and, in particular,

D@i;j ¼

Z Cj

ðiÞ @ðjÞ nij ds  /k wl ~

Z T i;j

@ðjÞ r/ðiÞ dx k wl

and

Q@i;j ¼

Z T i;j

@ðjÞ

ðiÞ

wk r/l dx 

Z Cj

@ðjÞ

ðiÞ

wk /l

ri;j~ nj ds:

P

~> ~ ~ j xj T xj

P 0.

h

79

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

 > ~ can be written as A ~ ¼ A þ B where Then it is still true that D@i;j ¼  Q@i;j and the complete system A

> X  @ ^ nþ1 Q@i;j M 1 j Qi;j pi

B : h Dt

j2Si \BðXÞ

"

> X  @ Q@i;j M 1 j Qi;j þ

~ : h Dt A

j2Si \BðXÞ

þ h Dt

#

X

Q>i;j M 1 j Qi;j

^ nþ1 p i

j2Si BðXÞ

X

^ nþ1 Q>i;j M 1 j Q}ði;jÞ;j p}ði;jÞ :

j2Si BðXÞ

It is easy to check that B is symmetric and at least positive semi-definite. We have to introduce now some types of boundary conditions in order to show that, if the pressure is specified on the ~ is positive definite. Let us rewrite x> Bx by including the external contribution and in boundary, the complete system A the form of the Eq. (45), namely

x> Bx ¼

Nj  h i X T @‘ðjÞ;j x‘ðjÞ þ T @ x j¼1

ext;j

>  h i T @‘ðjÞ;j x‘ðjÞ þ T @ x

ext;j

 ;

  1 where T @i;j ¼ Mj 2 Q@i;j and T @ x ext;j is a known external contribution that depends on the boundary conditions. In particular, if   the pressure is specified at the boundary, then T @ext;j ¼ T @‘ðjÞ;j and T @ x ext;j is a known quantity that in general is part of the   known right hand side vector. Since the external pressure is specified, then T @‘ðjÞ;j x‘ðjÞ þ T @ x ext;j ¼ 0 () x‘ðjÞ  xext;j . We take now x> Bx ¼ 0 that implicitly fixes xext ¼ 0. In this way x‘ðjÞ ¼ 0 8j 2 BðXÞ. Using the same reasoning on the matrix A we ~ is positive definite in this case. A possible way to specify the velocity at the boundary can conclude that x  0, and hence A is to neglect the jump contribution for the pressure at the boundary or equivalent, taken xext;j ¼ x‘ðjÞ 8j 2 BðXÞ. It is easy to ~ ¼ 0 for every x constant, and then the matrix A ~ is only check that if we have only this type of boundary conditions then x> Ax positive semi-definite. 2.8. Remarks on the stability In this section we discuss the stability of the proposed method. In particular we will show that if the discretization of the 2

d nonlinear convective-viscous term is stable, then the discrete solution of the velocity is L2 stable. Let v :¼ h ðuh ; v h Þ 2 ðV h Þ and m ph 2 V h be the numerical solution for velocity and pressure. Then the L2 stability property for the velocity reads

Z

d dt

1 2 u þ v 2h dx 6 0: 2 h

X

ð46Þ

By construction, the numerical solution satisfies the weak formulation of the continuity and momentum equation, namely:

Z

ðjÞ

Rj

þ

wk

Z

@ vh dx þ @t

Z @Rj

ðjÞ

T ‘ðjÞ;j

ðjÞ

n ds  wk Gh  ~

wk rph;‘ðjÞ dx þ

Z

Z Rj

rwðjÞ k  Fðv h ; rv h Þ dx

ðjÞ

T rðjÞ;j

wk rph;rðjÞ dx þ

Z Cj

  ðjÞ wk ph;rðjÞ  ph;‘ðjÞ ~ nstd ds ¼ 0

ð47Þ

and

" X Z j2Si

Cj

ðiÞ /k

vh  ~ne;j ds 

#

Z

ðiÞ /k

r

T i;j

 v h dx ¼ 0

ð48Þ

ðiÞ

for every i ¼ 1 . . . N i and j ¼ 1 . . . N j . Since f/k gk forms a basis for the polynomial space V m h jT i , we can take ph as a test function in (48). Denoting with ph;i the discrete solution of the pressure inside element number i we directly obtain from (48)

" X Z j2Si

ph;i v h  ~ ne;j ds 

Cj

#

Z T i;j

rph;i  v h dx ¼ 0:

ð49Þ

We can now sum over i ¼ 1 . . . N i in Eq. (49) and find

" Ni X Z X i¼1 j2Si

"Z Nj X j¼1

ph;i v h  ~ ne;j ds 

Cj

T ‘ðjÞ;j

rph;‘ðjÞ  v h dx þ

#

Z T i;j

rph;i  v h dx ¼ 0;

Z T rðjÞ;j

# Z   ~ rph;rðjÞ  v h dx þ ph;rðjÞ  ph;‘ðjÞ v h  nstd ds ¼ 0: Cj

ð50Þ

80

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

Using the same reasoning we can take as test function wk ¼ uh 2 V dh in the first component of Eq. (47) and wk ¼ v h 2 V dh in the second one. By summing the resulting two expressions we thus obtain

Z

vh 

Rj

þ

Z

@ v h dx þ @t

Z

vh  ðGh  ~nÞ ds  @Rj

v h  rph;‘ðjÞ dx þ

rv h : Fðv h ; rv h Þ dx

Rj

Z

T ‘ðjÞ;j

Z

v h  rph;rðjÞ dx þ

Z 

T rðjÞ;j

ph;rðjÞ  ph;‘ðjÞ



v h  ~nstd ds ¼ 0;

ð51Þ

Cj

where c ¼ A : B denotes the tensor operator c ¼ Aij Bij , written at the aid of the Einstein summation convention. Using Eq. (51) we obtain:

d dt

Z X

1 2 u þ v 2h dx ¼ 2 h

Z 

vh 

X

" X Z  ¼

  Nj Z  X @ @ v h dx ¼ v h  vh dx @t @t j¼1 Rj

Nj

j¼1



v h  ðGh  ~nÞ ds þ

#

rv h : Fðv h ; rv h Þ dx

Rj

@Rj

"Z Nj X j¼1

Z

v h  rph;‘ðjÞ dx þ

T ‘ðjÞ;j

Z T rðjÞ;j

v h  rph;rðjÞ dx þ

Z 

ph;rðjÞ  ph;‘ðjÞ



#

v h  ~nstd ds ;

ð52Þ

Cj

where the first sum on the right hand side contains the contribution of the nonlinear convective and viscous terms and the second sum, which contains the contribution of the pressure, vanishes thanks to (50). As a final result, we therefore obtain:

d dt

Z X

"Z # Z Nj X 1 2 2 u þ v h dx ¼  v h  ðGh  ~nÞ ds  rvh : Fðv h ; rv h Þ dx : 2 h @Rj Rj j¼1

ð53Þ

As a consequence of Eq. (53), the stability of our staggered DG method is determined by the stability of the discretization used for the nonlinear convective and the viscous term. We note that on the right hand side of Eq. (53) we have a standard DG discretization of these terms on the dual mesh, for which the standard DG stability results apply. In the literature there are stability results that prove that if the numerical flux Gh is an entropy flux, then the right hand side is actually less or equal zero for the convective part. More details on the stability of high order DG schemes for convective problems can be found in [56] for the scalar case and in [2] for first order systems. Furthermore, if we have no convective-viscous contributions, Eq. (53) states that the total kinetic energy of the system is exactly conserved. 3. Numerical test problems 3.1. Convergence test We consider a smooth steady state problem in order to measure the order of accuracy of the proposed method. For this purpose, the Navier–Stokes equations are first rewritten in cylindrical coordinates (r and u), with r 2 ¼ x2 þ y2 , tan u ¼ x=y, the radial velocity component ur and the angular velocity component uu . In order to derive an analytical solution we suppose a steady vortex-type flow with angular symmetry, i.e. @=@t ¼ 0; @=@ u ¼ 0 and ur ¼ 0. With these assumptions, the continuity equation is automatically satisfied and the system of incompressible Navier–Stokes equations reduces to

8 < @p ¼ u2u ; r

@r

: r @ 2 uu þ @ 2 uu  uu ¼ 0: @r r @r 2

ð54Þ

One can now recognize in the second equation of (54) a classical second order Cauchy Euler equation and so obtain two solutions for uu , namely:

uu ¼ c1 r; uu ¼

c1 r

ð55Þ ð56Þ

for every c1 2 R. The corresponding pressures read

p¼

c21 r 2 þ c2 ; 2

p ¼ 2

c21 þ c2 ; r2

ð57Þ

ð58Þ

respectively. In this section we set the boundary conditions in order to obtainp the non-trivial solution (56) and (58). Due to ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi the singularity of uu for r ¼ 0, let X ¼ Cð5Þ  Cð1Þ where CðrÞ ¼ fðx; yÞ 2 R2 j x2 þ y2 6 rg. As initial condition we impose

81

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

Eqs. (56)–(58) with c1 ¼ uu ð1Þ ¼ 2 and c2 ¼ 0. The exact velocity is imposed at the internal boundary while the exact pressure is specified at the external circle. The proposed algorithm is validated for several polynomial degrees p using successively refined grids. The chosen parameters for the numerical simulations are t end ¼ 0:75; h ¼ 1; m ¼ 105 ; the time step Dt is taken according to the CFL time restriction for the explicit discretization of the nonlinear convective term (35). The L2 error between the analytical and the numerical solution is computed as

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z ðpÞ ¼ ðph  pe Þ2 dx;

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z ðv Þ ¼ ðv h  v e Þ2 dx

X

ð59Þ

X

for the pressure and for the velocity vector field, respectively, where the subscript h indicates the numerical solution and e denotes the exact solution. Tables 1 and 2 show the L2 convergence rates for successive refinements of the grid, where OðpÞ and Oðv Þ represent the order of accuracy achieved for the pressure and the velocity field, respectively. The optimal convergence is reached up to p ¼ 2 while for p ¼ 3 the observable order of accuracy for the velocity vector field is closer to p þ 12 rather then p þ 1. 3.2. Womersley profiles In this section the proposed algorithm is verified against the exact solution for an oscillating flow in a rigid tube of length L. The unsteady flow is driven by a sinusoidal pressure gradient on the boundaries

~ pout ðtÞ  pinlet ðtÞ p ¼ eixt ; L q

ð60Þ

~ is the amplitude of the pressure gradient; q is the fluid density; x is the frequency of the oscillation; i indicates the where p imaginary unit; pinlet and pout are the inlet and outlet pressures, respectively. The analytical solution was derived by Womersley [82]. According to [82,41] no convective contribution is considered. By imposing Eq. (60) at the tube ends, the resulting unsteady velocity field is uniform in the axial direction and is given by

 3 3 2 J 0 afi2 ~ 1 p 41    5eixt ; ue ðx; y; tÞ ¼ 3 q ix J ai2

v e ðx; y; tÞ ¼ 0;

ð61Þ

0

pffiffiffi where f ¼ 2y=D is the dimensionless radial coordinate; D is the diameter of the tube; a ¼ D2 xm is a constant; and J 0 is the ~¼ zero-th order Bessel function of the first kind. For the present test we take X ¼ ½0:5; 1 ½0:2; 0:2; p 1000; q ¼ 1000; x ¼ 2p; h ¼ 0:6; and m ¼ 8:94 104 . The computational domain X is covered with a total number of N i ¼ 98 triangles and the time step size is chosen as Dt ¼ 0:01. The numerical results for p ¼ 3 are shown in Fig. 3 for several times at x ¼ 0:1. A good agreement between exact and numerical solution can be observed. 3.3. Blasius boundary layer Another classical test problem concerns the Blasius boundary layer. For the particular case of laminar stationary flow over a flat plate, a solution of Prandtl’s boundary layer equations was found by Blasius [10] and is determined by the solution of a third-order non-linear ODE, namely:

(

000

00

f þ ff ¼ 0; f ð0Þ ¼ 0;

0

f ð0Þ ¼ 0;

ð62Þ

0

lim f ðnÞ ¼ 1;

n!1

qffiffiffiffiffiffi 0 where n ¼ y 2um1x is the Blasius coordinate; f ¼ uu1 ; and u1 is the farfield velocity. The reference solution is computed here using a tenth-order DG ODE solver, see e.g. [38], together with a classical shooting method. In order to obtain the Blasius velocity profile in our simulations we consider a steady flow over a wedge-shaped object. As a result of the viscosity, a boundary layer appears along the obstacle. For the present test, we consider X ¼ ½0; 1 ½0:25; 0:25 and a wedge shape object with upper edge corresponding to the segment x ¼ ½0; 1. An initially uniform flow uðx; y; 0Þ ¼ u1 ¼ 1 , v ðx; y; 0Þ ¼ 0 and pðx; y; 0Þ ¼ 1 is imposed as initial condition, while an inflow boundary is imposed on the left and outflow boundary Table 1 Numerical convergence results for p ¼ 0 and p ¼ 1. Ni

p¼0

ðpÞ 124 496 1984 7936 31,744

7.902E01 5.026E01 2.982E01 1.659E01 8.797E02

p¼1

ð~ vÞ 1.095E00 7.086E01 4.502E01 2.797E01 1.714E01

OðpÞ

Oð~ vÞ

ðpÞ

– 0.7 0.8 0.8 0.9

– 0.6 0.7 0.7 0.7

3.944E01 8.830E02 2.325E02 6.207E03 1.615E03

ð~ vÞ 4.311E01 1.221E01 3.299E02 8.725E03 2.318E03

OðpÞ

Oð~ vÞ

– 2.2 1.9 1.9 1.9

– 1.8 1.9 1.9 1.9

82

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

Table 2 Numerical convergence results for p ¼ 2 and p ¼ 3. Ni

p¼2

p¼3

ðpÞ 124 496 1984 7936

9.366E02 1.054E02 1.193E03 1.438E04

ð~ vÞ 1.990E01 3.069E02 3.686E03 4.425E04

OðpÞ

Oð~ vÞ

ðpÞ

– 3.2 3.1 3.1

– 2.7 3.1 3.1

4.346E02 2.966E03 1.783E04 1.313E05

ð~ vÞ 9.317E02 8.027E03 7.153E04 5.997E05

OðpÞ

Oð~ vÞ

– 3.9 4.1 3.8

– 3.5 3.5 3.6

Exact solution Numerical solution (p=3)

0.2

y

0.1

0

-0.1

-0.2 -0.2

-0.1

0

0.1

u

0.2

Fig. 3. Comparison between the exact and the numerical solution for the Womersley profiles at times t ¼ 1:7; t ¼ 1:9; t ¼ 2:0; t ¼ 1:2, respectively, from left to right.

conditions are imposed on the other edges of the external box. Finally, no-slip wall boundary conditions are considered over the wedge shape object. We cover X with a total amount of N i ¼ 278 triangles and use h ¼ 1 and p ¼ 3. The resulting Blasius velocity profile is shown in Fig. 4, while the profile with respect to the Blasius coordinate n is shown in Fig. 5 in order to verify whether the obtained solution is self-similar with respect to n. A comparison between the numerical results presented here and the reference solution is depicted in Fig. 6 for x ¼ 0:4 and x ¼ 0:6. A good agreement between the reference solution and the numerical results obtained with the staggered semi-implicit DG scheme is obtained, despite the use of a very coarse grid. Note that the solution in terms of the Blasius coordinate n is

y

0.2

0

-0.2 0

x

0.5

1

Fig. 4. Computational domain used for the simulation of the Blasius boundary layer. The colors represent the horizontal velocity u.

83

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92 20

xi

15

10

5

0

0.2

0.4

x

0.6

0.8

1

Fig. 5. Velocity profile with respect to the Blasius coordinate n.

independent from x. The numerical solution is also verified to maintain the self-similar Blasius profile in the ðx; nÞ plane, see Fig. 5. 3.4. Lid-driven cavity flow We consider here another classical benchmark problem for the incompressible Navier–Stokes equations, namely the liddriven cavity problem. This test problem is solved numerically with the new staggered DG scheme on very coarse grids using a polynomial degree of p ¼ 3. Let X ¼ ½0:5; 0:5 ½0:5; 0:5, set velocity boundary conditions u ¼ 1 and v ¼ 0 on the top boundary (i.e. y ¼ 0:5) and impose no-slip wall boundary conditions on the other edges. As initial condition we take uðx; y; 0Þ ¼ v ðx; y; 0Þ ¼ 0. We use a grid with N i ¼ 73 triangles for Re ¼ 100; 400; 1000 and N i ¼ 359 triangles for Re ¼ 3200. A sketch of the main and dual grid is shown in Fig. 7. For the present test h ¼ 1; Dt is taken according to condition (35); and t end ¼ 150. According to [48], primary and corner vortices appear from Re ¼ 100 to Re ¼ 3200, a comparison of the velocities against the data presented in [48], as well as the streamline plots are shown in Fig. 8. A very good agreement is obtained in all cases, even if a very coarse grid has been used. 3.5. Backward-facing step. In this section, the numerical solution for the fluid flow over a backward-facing step is considered. For this test problem, both experimental and numerical results are available at several Reynolds numbers (see e.g. [1,40]). The computational domain X and the main notation are reported in Fig. 9. Reference u (x=0.4) u (x=0.6)

1.1 1 0.9 0.8 0.7

u

0.6 0.5 0.4 0.3 0.2 0.1 0 0

2

4

xi

6

8

Fig. 6. Numerical and reference solution for the Blasius boundary layer at x ¼ 0:4 and x ¼ 0:6.

84

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

0.4

y

0.2

0

-0.2

-0.4 -0.4

-0.2

0

x

0.2

0.4

Fig. 7. Main and dual grid used for the lid-driven cavity problem for Re ¼ 100; 400; 1000.

The fluid flow is driven by a pressure gradient imposed at the left and the right ends of the computational domain. On all the other boundaries, no-slip wall boundary conditions are imposed. According to [1], we take Re ¼ DU m where D ¼ 2hin ; U is the mean inlet velocity; m is the kinematic viscosity. The computational domain is covered with a total number of N i ¼ 260 triangles with characteristic size h ¼ 0:2 for x 6 5 and h ¼ 0:48 for x > 5 (see Fig. 9). Finally we use p ¼ 3; h ¼ 1 and Dt is the one given by the CFL condition for the nonlinear convective term; tend ¼ 80 s. Fig. 11 shows the vortices generated at different Reynolds numbers, while in Fig. 10 the main recirculation point X1 is compared with experimental data given by Armaly [1], and the explicit second-order upwind finite difference scheme introduced in [9]. A good agreement with the experimental data is shown up to Re ¼ 316 but, according to [1], the experiment becomes three dimensional for Re > 400, so the comparison can be done only up to Re ¼ 400. Indeed, one can see in Fig. 11 how the secondary vortex occurs for Re ¼ 426, while in the experiments it appears at higher Reynolds numbers (see e.g. [1]). 3.6. Rotational flow past a circular half-cylinder Here we consider a rotational flow past a circular half-cylinder. A comparison between numerical and exact analytical solution is possible for incompressible and inviscid fluid, i.e. here we set m ¼ 0. We use the computational setup of Feistauer pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and Kucera [42], hence X ¼ ½5; 5 ½0; 5  f x2 þ y2 6 0:5g; as boundary conditions we impose the velocity at the left boundary; homogeneous Neumann boundary conditions on the top and right boundaries and inviscid wall at the bottom and the surface of the half-cylinder. The farfield velocity field is given by u ¼ y and v ¼ 0. The exact analytical solution to this problem was found by Fraenkel in [45]. For the present test we choose p ¼ 3; Dt is set according to (35) and we cover X with Ni ¼ 800 triangles, using only 6 triangles to describe the half-cylinder. Curved isoparametric elements are considered in order to represent the geometry of the half-cylinder properly. As initial conditions we impose pðx; y; 0Þ ¼ 1; uðx; y; 0Þ ¼ y and v ðx; y; 0Þ ¼ 0. Two vortices appear near the half-cylinder (see Fig. 12 left), while a comparison between analytical and numerical velocity magnitude on the cylinder surface (i.e. r ¼ 0:5) is shown on the right of Fig. 12. A good agreement between analytical and numerical results is obtained also with a very coarse grid. An important remark is that for this test problem the use of isoparametric elements is crucial, as previously shown for inviscid flow past a circular cylinder by Bassi and Rebay [5]. 3.7. Flow over a circular cylinder In this section we consider the flow over a circular cylinder. Also in this case, the use of the isoparametric approach is mandatory to represent the geometry of the cylinder wall, see [5,71]. In particular, two cases are considered: first, an inviscid flow around the cylinder is assumed in order to obtain a steady potential flow; finally, the complete viscous case is considered in order topget the unsteady von Karman vortex street. For the first case a sufficiently large domain ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X ¼ ½8; 8 ½8; 8  f x2 þ y2 6 1g is employed. The exact solution for this case is known and reads:

 1 ur ðr; uÞ ¼ u

! R2c cosðuÞ; r2

 1þ uu ðr; uÞ ¼ u

! R2c sinðuÞ; r2

85

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92 vRef(Ghia) uRef(Ghia) u(p=3) v(p=3)

Re100

1.2

Re100

1

0.4

Velocities

0.8

0.2

0.6

0

y

0.4 0.2

-0.2

0 -0.2 -0.4

-0.4 -0.4

-0.2

0

xi

0.2

0.4

Re400

-0.4

-0.2

0

x

0.2

0.4

v RefRe400

uRefRe400 u(p=3) u(p=3)

1.2

Re400

1 0.8

0.4

Velocities

0.6 0.2

0.4

y

0.2

0

0 -0.2

-0.2 -0.4 -0.6

-0.4 -0.4

-0.2

0

xi

0.2

0.4

Re1000

-0.4

-0.2

v RefRe1000 uRefRe1000 v(p=3) u(p=3)

1.2

0

x

0.2

0.4

0.2

0.4

Re1000

1

0.4

0.8

Velocities

0.6

0.2

y

0.4 0.2 0

-0.2

-0.2 -0.4 -0.6

0

-0.4 -0.4

-0.2

0

xi

0.2

0.4

-0.2

v Ref(Ghia) uRef(Ghia) u(p=3) v(p=3)

Re3200

1.2

-0.4

0

x

Re3200

1

0.4

0.6

0.2

0.4 0.2

y

Velocities

0.8

0

-0.2

25

-0.4 -0.6 -0.8

0 -0.2 -0.4

-0.4

-0.2

0

xi

0.2

0.4

-0.4

-0.2

0

x

0.2

0.4

Fig. 8. Velocity profiles (left) and streamlines (right) at several Reynolds numbers for the lid-driven cavity problem.

86

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

h_in

Second Recirculation

h Main Recirculation

X1

Fig. 9. Grid and main notation used for the backward-facing step problem.

14

Experimental data Staggered semi-implicit DG scheme (p=3) Numerical solution by Erturk

12

10

X1

8

6

4

2

0

0

100

200

300

Re

400

500

600

700

Fig. 10. Comparison of the experimental data of Armaly et al. [1] with the numerical results obtained with the present semi-implicit staggered DG scheme and the numerical solution obtained in [40] for the reattachment point X1 in the backward-facing step problem.

! 1 2 2R2c R4c  p¼ u cosð2uÞ  4 ; 2 r2 r

ð63Þ

 is the inflow velocity; Rc is the cylinder radius; ur and uu are the radial and angular components of the velocity, where u  ; 0Þ is used, while the exact velocity distribution is taken as the external respectively. An initial condition v ðx; y; 0Þ ¼ ðu boundary condition. An inviscid wall boundary condition is imposed on the cylinder. For the present test  ¼ 0:01; Rc ¼ 1; m ¼ 0; p ¼ 3; h ¼ 0:6; Dt is the one taken according to the CFL restriction (35); t end ¼ 10. The domain X u is covered with a total number of N i ¼ 1464 triangles and an isoparametric approach is considered to represent the cylinder wall properly. Fig. 13 shows the streamlines and the pressure contours obtained at t ¼ 10 as well as the comparison between exact and numerical solution at several radii. A very good agreement between exact and numerical solution is observed. We consider now the fully viscous case in order to show the formation of the von Karman vortex street. Two domains are considered here: X1 ¼ ½20; 80 ½20; 20 covered with a N i ¼ 1702 triangles; and X2 ¼ ½5; 30 ½10; 10 covered with a  ; 0Þ; and u  ¼ 0:5. Different viscosity coefficients are used in order N i ¼ 1706 triangles. As initial condition we set v ðx; y; 0Þ ¼ ðu  ; 0Þ is to obtain different Reynolds numbers. For the present test we use Dt according to (35); p ¼ 3; h ¼ 0:6. The velocity ðu prescribed at the left boundary while homogeneous Neumann boundary conditions are imposed on the other external edge of the domains. Finally viscous wall boundary condition is imposed on the cylinder surface. Fig. 14 shows the obtained relationship between the Strouhal number, computed as St ¼ 2rf , the numerical results given u1 by Qu et al. (see [64]) and the experimental law given in [81]. The simulations are performed on the domain X1 . The numerical results fit well the experimental data and the numerical reference solution up to Re ¼ 150. Better results can be obtained by further enlarging the computational domain. The velocity field and the vorticity show different structures when low and high Reynolds numbers are considered. The vorticity contours are shown in Figure 15 for Re ¼ 50 and Re ¼ 125 at time t ¼ 500. In the case of Re ¼ 125 the von Karman vortex street is fully developed while, for Re ¼ 50, the two initial vortices remain present behind the cylinder for a longer

87

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

0

1

2

3

4

5

0

1

2

3

4

5

0

1

2

3

4

5

0

1

2

3

4

5

0

1

2

3

4

5

0

1

2

3

4

5

x

x

x

x

x

x

6

7

8

9

10

6

7

8

9

10

6

7

8

9

10

6

7

8

9

10

6

7

8

9

10

6

7

8

9

10

Fig. 11. Streamlines at Reynolds numbers Re ¼ 44; 113; 250; 316; 426 and 633 from top to bottom.

time. This is due to the low value of the Reynolds number, taken close to the limit of Re ¼ 40 for the generation of the vortex street. The time evolution of the generation of the von Karman vortex street is presented at several times for Re ¼ 200 on X2 in Fig. 16. Finally, in Fig. 17 we report a comparison between the computational time needed per time step for the main parts of the algorithm presented in this paper up to the time t ¼ 10 s using Re ¼ 100 on X1 if we employ a GMRES method or the cheaper

88

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92 Numerical Exact solution

1 1.5

0.8

0.6

y

Velocity

1

0.5

0.4

0.2 0 0 -0.5 -2

-1.5

-1

-0.5

0

0.5

x

1

1.5

2

-0.2 0

0.5236 1.0472 1.5708 2.0944 2.618 3.1416

φ

Fig. 12. Rotational inviscid flow past a circular half-cylinder. Left: Streamlines. Right: Analytical and numerical results for r ¼ 0:5. The vertical lines show the dimension of the six curved elements that cover the half-cylinder.

p 3

1.00005

1.00004 1.00003 1.00002 1.00001

2

Numerical Reference r=1.0 Reference r=1.5 Reference r=2.0

0.025

1

y

0

-1

-2

-3 -2

0

2

x

Numerical Reference r=1.0 Reference r=1.5 Reference r=2.0

0.015 0.01

0.015 0.01 0.005 0 -0.005

-2

-1

0

Theta

1

2

3

Numerical Reference r=1.0 Reference r=1.5 Reference r=2.0

1.0001

p

1

0

0.99995

-0.005

0.9999

-0.01

0.99985

-0.015

-3

1.00005

0.005

v

0.02

u

0.99999 0.99998 0.99997 0.99996 0.99995 0.99994 0.99993 0.99992 0.99991 0.9999 0.99989 0.99988 0.99987 0.99986 0.99985

1

-3

-2

-1

0

Theta

1

2

3

0.9998

-3

-2

-1

0

Theta

1

2

3

Fig. 13. Steady flow of an inviscid incompressible fluid around a circular cylinder. On the top left: streamlines and pressure contours at tend ¼ 10; Numerical and exact solution at r ¼ 1:0; r ¼ 1:5 and r ¼ 2:0 for the velocity components u; v and pressure p from top right to the bottom right, respectively.

CG method for the solution of the linear system. Note that since our particular semi-implicit DG discretization of the incompressible Navier–Stokes equations on staggered grids leads to a symmetric and positive-definite linear system, we can employ the CG method. This is not always the case for DG schemes applied to the incompressible Navier–Stokes equations since some formulations may also lead to non-symmetric linear systems.

89

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92 Williamson and Brown (Experimental) Datas from Lixia Qu et al Staggered DG p=3

0.2

St

0.18

0.16

0.14

0.12 60

80

100

120

140

160

180

Re Fig. 14. Strouhal–Reynolds number relationship for the present method, the method of Qu et al. [64] and experimental data of Williamson and Brown [81].

0.178378 0.151351 0.124324 0.0972973 0.0702703 0.0432432 0.0162162 -0.0108108 -0.0378378 -0.0648649 -0.0918919 -0.118919 -0.145946 -0.172973 -0.2

0

20

x

40

60

0.445946 0.378378 0.310811 0.243243 0.175676 0.108108 0.0405405 -0.027027 -0.0945946 -0.162162 -0.22973 -0.297297 -0.364865 -0.432432 -0.5

0

5

10

15

20

25

30

35 x

40

45

50

55

60

65

70

75

Fig. 15. Dual mesh and vorticity contours of the von Karman vortex street generated at time t ¼ 500 for Re ¼ 50 (top) and Re ¼ 125 (bottom).

The time required to compute the convective-viscous term represents in the second case the main computational effort. Using the GMRES algorithm the computational time needed to solve the linear system increases a lot compared to the CG method and becomes the main cost of the algorithm. In particular, the mean time to solve the system using the GMRES algorithm is, for this test, 6.2 s while using the CG method is only about 1.0 s. For all tests, the tolerance for solving the linear system was set to tol ¼ 1012 . We underline that for a fair comparison of the two methods, no preconditioners have been used and that faster convergence can be obtained by using a proper preconditioner for each iterative solver. 4. Conclusions A new, spatially high order accurate semi-implicit DG scheme for the solution of the incompressible Navier–Stokes equations on staggered unstructured curved meshes has been proposed. The high order of accuracy in space was verified and compared with reference solutions for polynomial degrees up to p ¼ 3. The numerical results agree very well with

90

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

10

10

x

x

20

30

10

20

30

10

x

x

20

30

20

30

Fig. 16. Temporal evolution of the vorticity profile for t ¼ 15; t ¼ 30; t ¼ 50; t ¼ 75 from top left to bottom right at Re ¼ 200.

CGM Total Time CGM Fv Time CGM time GMRES Total Time GMRES Fv Time GMRES time

20

CGM GMRES

3500

3000 15

Iterations

time

2500

10

2000

1500 5 1000

0

200

400

n

600

800

500

200

400

n

600

800

Fig. 17. Left: comparison between the CPU time for the GMRES method and the CG method. Right: total number of iterations for the GMRES against the number of iterations for the CG method.

the reference data for all test cases considered in this paper. The proposed numerical method reduces to a classical semi-implicit finite-volume and finite-difference scheme on staggered unstructured meshes [21] for p ¼ 0 if the dual grid is chosen orthogonal to the primary mesh, i.e. if the dual mesh is based on the circumcenters of the primary triangles, and if the pressure jump term in the momentum equation is multiplied with a factor of 12. The same comment also applies to the case p ¼ 0 for the method presented in [71] for the shallow water equations. Furthermore, the use of matrices that depend only on the geometry and on the polynomial degree and hence can be precomputed before runtime, leads to a computationally efficient scheme. In addition, the resulting main matrix results symmetric and positive definite for appropriate boundary conditions. This allows to use fast iterative methods for the solution of the sparse linear system with a significant gain in terms of computational time. Future research will concern the extension of the scheme to high order of accuracy also in time using a space–time DG approach as well as the extension to the fully three-dimensional case on unstructured tetrahedral meshes. Acknowledgments The research presented in this paper was funded by the European Research Council (ERC) under the European Union’s Seventh Framework Programme (FP7/2007–2013) within the research project STiMulUs, ERC Grant agreement No. 278267.

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

91

References [1] B.F. Armaly, F. Durst, J.C.F. Pereira, B. Schonung, Experimental and theoretical investigation on backward-facing step flow, J. Fluid Mech. 127 (1983) 473–496. [2] T. Barth, Pierre Charrier, Energy stable flux formulas for the discontinuous Galerkin discretization of first-order nonlinear conservation laws, Technical Report NAS-01-001, NASA, 2001. [3] F. Bassi, A. Crivellini, D.A. Di Pietro, S. Rebay, On a robust discontinuous Galerkin technique for the solution of compressible flow, J. Comput. Phys. 218 (2006) 208–221. [4] F. Bassi, A. Crivellini, D.A. Di Pietro, S. Rebay, An implicit high-order discontinuous Galerkin method for steady and unsteady incompressible flows, Comput. Fluids 36 (2007) 1529–1546. [5] F. Bassi, S. Rebay, A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations, J. Comput. Phys. 131 (1997) 267–279. [6] M. Ben-Artzi, J. Falcovitz, A second-order Godunov-type scheme for compressible fluid dynamics, J. Comput. Phys. 55 (1984) 1–32. [7] A. Bermudez, A. Dervieux, J.A. Desideri, M.E. Vazquez, Upwind schemes for the two-dimensional shallow water equations with variable depth using unstructured meshes, Comput. Methods Appl. Mech. Eng. 155 (1998) 49–72. [8] A. Bermúdez, J.L. Ferrín, L. Saavedra, M.E. Vázquez-Cendón, A projection hybrid finite volume/element method for low-Mach number flows, J. Comput. Phys. 271 (2014) 360–378. [9] F. Biagio, Calculation of laminar flows with second-order schemes and collocated variable arrangement, Int. J. Numer. Methods Fluids 26 (1998) 887– 905. [10] H. Blasius, Grenzschichten in Flüssigkeiten mit kleiner Reibung, Z. Math. Phys. 56 (1908) 1–37. [11] A.N. Brooks, T.J.R. Hughes, Stream-line upwind/Petrov Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equation, Comput. Methods Appl. Mech. Eng. 32 (1982) 199–259. [12] L. Brugnano, V. Casulli, Iterative solution of piecewise linear systems, SIAM J. Sci. Comput. 30 (2007) 463–472. [13] L. Brugnano, V. Casulli, Iterative solution of piecewise linear systems and applications to flows in porous media, SIAM J. Sci. Comput. 31 (2009) 1858– 1873. [14] L. Brugnano, A. Sestini, Iterative solution of piecewise linear systems for the numerical solution of obstacle problems, J. Numer. Anal. Ind. Appl. Math. 6 (2012) 67–82. [15] V. Casulli, A semi-implicit finite difference method for non-hydrostatic free-surface flows, Int. J. Numer. Methods Fluids 30 (1999) 425–440. [16] V. Casulli, A high-resolution wetting and drying algorithm for free-surface hydrodynamics, Int. J. Numer. Methods Fluids 60 (2009) 391–408. [17] V. Casulli, A semi-implicit numerical method for the free-surface Navier–Stokes equations, Int. J. Numer. Methods Fluids 74 (2014) 605–622. [18] V. Casulli, E. Cattani, Stability, accuracy and efficiency of a semi-implicit method for three-dimensional shallow water flow, Comput. Math. Appl. 27 (1994) 99–112. [19] V. Casulli, R.T. Cheng, Semi-implicit finite difference methods for three-dimensional shallow water flow, Int. J. Numer. Methods Fluids 15 (1992) 629– 648. [20] V. Casulli, G.S. Stelling, Semi-implicit subgrid modelling of three-dimensional free-surface flows, Int. J. Numer. Methods Fluids 67 (2011) 441–449. [21] V. Casulli, R.A. Walters, An unstructured grid, three-dimensional model based on the shallow water equations, Int. J. Numer. Methods Fluids 32 (2000) 331–348. [22] V. Casulli, P. Zanolli, High resolution methods for multidimensional advection–diffusion problems in free-surface hydrodynamics, Ocean Modell. 10 (2005) 137–151. [23] V. Casulli, P. Zanolli, Iterative solutions of mildly nonlinear systems, J. Comput. Appl. Math. 236 (2012) 3937–3947. [24] A.J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comput. Phys. 2 (1967) 12–26. [25] A.J. Chorin, Numerical solution of the Navier–Stokes equations, Math. Comput. 23 (1968) 341–354. [26] E.T. Chung, C.S. Lee, A staggered discontinuous Galerkin method for the convection–diffusion equation, J. Numer. Math. 20 (2012) 1–31. [27] E.T. Chung, P. Ciarlet, T.F. Yu, Convergence and superconvergence of staggered discontinuous Galerkin methods for the three-dimensional Maxwell’s equations on Cartesian grids, J. Comput. Phys. 235 (2013) 14–31. [28] E.T. Chung, B. Engquist, Optimal discontinuous Galerkin methods for wave propagation, SIAM J. Numer. Anal. 44 (2006) 2131–2158. [29] E.T. Chung, B. Engquist, Optimal discontinuous Galerkin methods for the acoustic wave equation in higher dimensions, SIAM J. Numer. Anal. 47 (2009) 3820–3848. [30] E.T. Chung, H.H. Kim, O.B. Widlund, Two-level overlapping Schwarz algorithms for a staggered discontinuous Galerkin method, SIAM J. Numer. Anal. 51 (2013) 47–67. [31] B. Cockburn, C.W. Shu, The local discontinuous Galerkin method for time-dependent convection diffusion systems, SIAM J. Numer. Anal. 35 (1998) 2440–2463. [32] B. Cockburn, C.W. Shu, The Runge–Kutta discontinuous Galerkin method for conservation laws. V: Multidimensional systems, J. Comput. Phys. 141 (1998) 199–224. [33] B. Cockburn, C.W. Shu, Runge–Kutta discontinuous Galerkin methods for convection-dominated problems, J. Sci. Comput. 16 (2001) 173–261. [34] A. Crivellini, V. D’Alessandro, F. Bassi, High-order discontinuous Galerkin solutions of three-dimensional incompressible RANS equations, Comput. Fluids 81 (2013) 122–133. [35] V. Dolejsi, Semi-implicit interior penalty discontinuous Galerkin methods for viscous compressible flows, Commun. Comput. Phys. 4 (2008) 231–274. [36] V. Dolejsi, M. Feistauer, A semi-implicit discontinuous Galerkin finite element method for the numerical solution of inviscid compressible flow, J. Comput. Phys. 198 (2004) 727–746. [37] V. Dolejsi, M. Feistauer, J. Hozman, Analysis of semi-implicit DGFEM for nonlinear convection–diffusion problems on nonconforming meshes, Comput. Methods Appl. Mech. Eng. 196 (2007) 2813–2827. [38] M. Dumbser, Arbitrary high order PNPM schemes on unstructured meshes for the compressible Navier–Stokes equations, Comput. Fluids 39 (2010) 60–76. [39] M. Dumbser, V. Casulli, A staggered semi-implicit spectral discontinuous Galerkin scheme for the shallow water equations, Appl. Math. Comput. 219 (15) (2013) 8057–8077. [40] E. Erturk, Numerical solutions of 2d steady incompressible flow over a backward-facing step. Part I: High reynolds number solutions, Comput. Fluids 37 (2008) 633–655. [41] F. Fambri, M. Dumbser, V. Casulli. An efficient semi-implicit method for three-dimensional non-hydrostatic flows in compliant arterial vessels, Int. J. Numer. Methods Biomed. Eng., in press, http://dx.doi.org/10.1002/cnm.2651. [42] M. Feistauer, V. Kucera, On a robust discontinuous Galerkin technique for the solution of compressible flow, J. Comput. Phys. 224 (2007) 208–221. [43] E. Ferrer, R.H.J. Willden, A high order discontinuous Galerkin finite element solver for the incompressible Navier–Stokes equations, Comput. Fluids 46 (2011) 224–230. [44] M. Fortin, Old and new finite elements for incompressible flows, Int. J. Numer. Methods Fluids 1 (1981) 347–364. [45] L. Fraenkel, On corner eddies in plane inviscid shear flow, J. Fluid Mech. 11 (1961) 400–406. [46] G. Gassner, F. Lörcher, C.D. Munz, A contribution to the construction of diffusion fluxes for finite volume and discontinuous Galerkin schemes, J. Comput. Phys. 224 (2007) 1049–1063. [47] F.X. Giraldo, M. Restelli, High-order semi-implicit time-integrators for a triangular discontinuous Galerkin oceanic shallow water model, Int. J. Numer. Methods Fluids 63 (2010) 1077–1102.

92

M. Tavelli, M. Dumbser / Applied Mathematics and Computation 248 (2014) 70–92

[48] U. Ghia, K.N. Ghia, C. T. Shin. High-Re solutions for incompressible flow using Navier-Stokes equations, multigrid method, J.Comput. Phys. (48) (1982) 387–411. [49] S. Gottlieb, C.W. Shu, Total variation diminishing Runge–Kutta schemes, Math. Comput. 67 (1998) 73–85. [50] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface, Phys. Fluids 8 (1965) 2182– 2189. [51] M.R. Hestenes, E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Natl. Bur. Stand. 49 (1952) 409–436. [52] J.G. Heywood, R. Rannacher, Finite element approximation of the nonstationary Navier–Stokes Problem. I: Regularity of solutions and second order error estimates for spatial discretization, SIAM J. Numer. Anal. 19 (1982) 275–311. [53] J.G. Heywood, R. Rannacher, Finite element approximation of the nonstationary Navier–Stokes Problem. III: Smoothing property and higher order error estimates for spatial discretization, SIAM J. Numer. Anal. 25 (1988) 489–512. [54] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [55] T.J.R. Hughes, M. Mallet, M. Mizukami, A new finite element formulation for computational fluid dynamics. II: Beyond SUPG, Comput. Methods Appl. Mech. Eng. 54 (1986) 341–355. [56] G. Jiang, C.W. Shu, On a cell entropy inequality for discontinuous Galerkin methods, Math. Comput. 62 (1994) 531–538. [57] H.H. Kim, E.T. Chung, C.S. Lee, A staggered discontinuous Galerkin method for the stokes system, SIAM J. Numer. Anal. 51 (2013) 3327–3350. [58] B. Klein, F. Kummer, M. Oberlack, A SIMPLE based discontinuous Galerkin solver for steady incompressible flows, J. Comput. Phys. 237 (2013) 235–250. [59] Y.J. Liu, C.W. Shu, E. Tadmor, M. Zhang, Central discontinuous Galerkin methods on overlapping cells with a non-oscillatory hierarchical reconstruction, SIAM J. Numer. Anal. 45 (2007) 2442–2467. [60] Y.J. Liu, C.W. Shu, E. Tadmor, M. Zhang, L2-stability analysis of the central discontinuous Galerkin method and a comparison between the central and regular discontinuous Galerkin methods, Math. Model. Numer. Anal. 42 (2008) 593–607. [61] 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. [62] V.S. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, 1980. [63] V.S. Patankar, B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, Int. J. Heat Mass Transfer 15 (1972) 1787–1806. [64] L. Qu, C. Norberg, L. Davidson, S.H. Peng, F. Wang, Quantitative numerical analysis of flow past a circular cylinder at reynolds number between 50 and 200, J. Fluids Struct. 39 (2013) 347–370. [65] S. Rhebergen, B. Cockburn, A spacetime hybridizable discontinuous Galerkin method for incompressible flows on deforming domains, J. Comput. Phys. 231 (2012) 4185–4204. [66] S. Rhebergen, B. Cockburn, Jaap J.W. van der Vegt, A space–time discontinuous Galerkin method for the incompressible Navier–Stokes equations, J. Comput. Phys. 233 (2013) 339–358. [67] V.V. Rusanov, Calculation of interaction of non-steady shock waves with obstacles, J. Comput. Math. Phys. USSR 1 (1961) 267–279. [68] Y. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869. [69] K. Shahbazi, P.F. Fischer, C.R. Ethier, A high-order discontinuous Galerkin method for the unsteady incompressible Navier–Stokes equations, J. Comput. Phys. 222 (2007) 391–407. [70] C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys. 77 (1988) 439–471. [71] M. Tavelli, M. Dumbser, A high order semi-implicit discontinuous Galerkin method for the two dimensional shallow water equations on staggered unstructured meshes, Appl. Math. Comput. 234 (2014) 623–644. [72] M. Tavelli, M. Dumbser, V. Casulli, High resolution methods for scalar transport problems in compliant systems of arteries, Appl. Numer. Math. 74 (2013) 62–82. [73] C. Taylor, P. Hood, A numerical solution of the Navier–Stokes equations using the finite element technique, Comput. Fluids 1 (1973) 73–100. [74] V.A. Titarev, E.F. Toro, ADER schemes for three-dimensional nonlinear hyperbolic systems, J. Comput. Phys. 204 (2005) 715–736. [75] E.F. Toro, A. Hidalgo, M. Dumbser, FORCE schemes on unstructured meshes. I: Conservative hyperbolic systems, J. Comput. Phys. 228 (2009) 3368– 3389. [76] E.F. Toro, V.A. Titarev, Solution of the generalized Riemann problem for advection–reaction equations, Proc. R. Soc. London (2002) 271–281. [77] G. Tumolo, L. Bonaventura, M. Restelli, A semi-implicit, semi-Lagrangian, p-adaptive discontinuous Galerkin method for the shallow water equations, J. Comput. Phys. 232 (2013) 46–67. [78] J. van Kan, A second-order accurate pressure correction method for viscous incompressible flow, SIAM J. Sci. Stat. Comput. 7 (1986) 870–891. [79] R. Verfürth, Finite element approximation of incompressible Navier–Stokes equations with slip boundary condition II, Numer. Math. 59 (1991) 615– 636. [80] R.A. Walters, V. Casulli, A robust finite element model for hydrostatic surface water flows, Commun. Numer. Methods Eng. 14 (1998) 931–940. pffiffiffiffiffi ffi [81] C.H.K. Williamson, G.L. Brown, A series in 1= Re to represent the Strouhal–Reynolds number relationship of the cylinder wake, J. Fluids Struct. 12 (1998) 1073–1085. [82] J. Womersley, Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known, J. Physiol. 127 (1955) 553–563.