A staggered space–time discontinuous Galerkin method for the incompressible Navier–Stokes equations on two-dimensional triangular meshes

Computers & Fluids 119 (2015) 235–249

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

A staggered space–time discontinuous Galerkin method for the incompressible Navier–Stokes equations on two-dimensional triangular meshes Maurizio Tavelli a, Michael Dumbser b,⇑ a b

Department of Mathematics, University of Trento, Via Sommarive 14, I-38050 Trento, Italy 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

Article history: Received 21 April 2015 Received in revised form 25 June 2015 Accepted 3 July 2015 Available online 8 July 2015 Keywords: Incompressible Navier–Stokes equations Semi-implicit space–time discontinuous Galerkin scheme Staggered unstructured meshes Space–time pressure correction method High order accuracy in space and time Isoparametric finite elements for curved boundaries

a b s t r a c t In this paper we propose a novel arbitrary high order accurate semi-implicit space–time discontinuous Galerkin method for the solution of the two dimensional incompressible Navier–Stokes equations on staggered unstructured triangular meshes. Isoparametric finite elements are used to take into account curved domain boundaries. The discrete pressure is defined on the primal triangular grid and the discrete velocity field is defined on an edge-based staggered dual grid. While staggered meshes are state of the art in classical finite difference approximations of the incompressible Navier–Stokes equations, their use in the context of high order DG schemes is novel and still quite rare. Formal substitution of the discrete momentum equation on the dual grid into the discrete continuity equation on the primary grid yields a very sparse four-point block system for the scalar pressure, which is conveniently solved with a matrix-free GMRES algorithm. Note that the same space–time DG scheme on a collocated grid would lead to ten non-zero blocks per element, since substituting the discrete velocity into the discrete continuity equation on a collocated mesh would involve also neighbors of neighbors. From numerical experiments we find that our linear system is well-behaved and that the GMRES method converges quickly even without the use of any preconditioner, which is a unique feature in the context of high order implicit DG schemes. A very simple and efficient Picard iteration is then used in order to derive a space–time pressure correction algorithm that achieves also high order of accuracy in time, which is in general a non-trivial task in the context of high order discretizations for the incompressible Navier–Stokes equations. The flexibility and accuracy of high order space–time DG methods on curved unstructured meshes allows to discretize even complex physical domains with very coarse grids in both, space and time. The use of a staggered grid allows to avoid the use of Riemann solvers in several terms of the discrete equations and significantly reduces the total stencil size of the linear system that needs to be solved for the pressure. The proposed method is verified for approximation polynomials of degree up to p ¼ 4 in space and time by solving a series of typical numerical test problems and by comparing the obtained numerical results with available exact analytical solutions or other numerical reference data. To the knowledge of the authors, this is the first time that a space–time discontinuous Galerkin finite element method is presented for the incompressible Navier–Stokes equations on staggered unstructured grids. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The discretization of the incompressible Navier–Stokes equations was mainly carried out in the past using finite difference methods [1–4] as well as continuous finite element schemes

⇑ Corresponding author. E-mail addresses: [email protected] (M. Tavelli), [email protected] (M. Dumbser). http://dx.doi.org/10.1016/j.compfluid.2015.07.003 0045-7930/Ó 2015 Elsevier Ltd. All rights reserved.

[5–11]. On the contrary, the construction of high order discontinuous Galerkin (DG) finite element methods for the incompressible Navier–Stokes equations is still a very active topic of ongoing research. Obtaining high order of accuracy also in time represents an important goal in order to achieve accurate results for unsteady problems. Several high order DG methods for the incompressible Navier– Stokes equations have been recently presented in the literature, see for example [12–19], without pretending completeness. An

236

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

alternative is the DG scheme proposed by Bassi et al. in [20], which is based on an extension of the technique of artificial compressibility that was originally introduced by Chorin in the finite difference context [21,22]. Another very well known approach to discretize general convection–diffusion equations in the context of hp discontinuous Galerkin finite element methods is the one proposed by Baumann and Oden in [23,24]. A unified analysis of several variants of the DG method applied to an elliptic model problem has been provided by Arnold et al. in [25]. We also would like to mention recent works on semi-implicit DG schemes, such as the ones presented in [26–30], to which our approach is indirectly related. While the use of staggered grids is a very common practice in the finite difference community, its use is not so widespread in the development of high order DG schemes. The first staggered DG schemes, based on a vertex-based dual grid, have been proposed in [31,32]. Other recent high order staggered DG schemes that use an edge-based dual grid have been forwarded in [33–35]. The advantage in using edge-based staggered grids is that they allow to improve significantly the sparsity pattern of the final linear system that has to be solved for the pressure. Very recently, a staggered semi-implicit DG scheme for the solution of the two dimensional shallow water equations was presented in [35,36] and then extended in [37] to the incompressible Navier–Stokes equations. The method presented in [37] is in principle of arbitrary high order of accuracy in space, while it reaches only second order in time. Consequently, it does not allow to recover high order accurate results for fully unsteady solutions. In this paper we propose a new method that is based on the general ideas put forward in [35–37], but which is also able to reach high order of accuracy in time. For this purpose we construct an arbitrary high order accurate staggered space–time discontinuous Galerkin finite element scheme. By relying on staggered grids we follow the classical philosophy of staggered finite difference schemes for the incompressible Navier–Stokes equations and for the free surface shallow water and Navier–Stokes equations, see [1–4,38–44]. In the context of staggered finite difference schemes we also would like to mention the so-called multiple pressure variables approach (MPV) [45–47], which is based on the asymptotic analysis of the compressible Navier–Stokes equations and is able to preserve also their incompressible limit. Our staggered semi-implicit space–time DG method proposed in this paper can be seen as a natural extension of the staggered semi-implicit DG scheme proposed in [37] to arbitrary high order of accuracy also in time. However, we emphasize that this extension is not straightforward for the complete convective-viscous problem. In the staggered DG scheme presented in [37], the discrete pressure is defined on the control volumes of the primal triangular mesh, while the discrete velocity vector field is defined on an edge-based, quadrilateral dual mesh. In the proposed staggered space–time DG scheme, the spatial control volumes are simply extended to the corresponding space–time control volumes by using the tensor product of the spatial control volume with the time interval of each time step, hence leading to triangular base prisms for the primal mesh and to quadrilateral base prisms for the dual mesh. For an overview of standard space–time DG schemes on collocated grids, see [16,17,48–50], and for a very early implementation of space–time finite element schemes in the context of moving meshes the reader is referred to [51]. The nonlinear convective terms are discretized explicitly by using a standard DG scheme [52–54] based on the local Lax-Friedrichs (Rusanov) flux [55], while the viscous terms are discretized implicitly using a fractional step method.1 The DG

1 Note that high order in time is obtained later by a Picard iteration, which embraces the entire scheme in each time step.

discretization of the viscous fluxes is based on the formulation of Gassner et al. [56], who obtained the viscous numerical flux from the solution of the Generalized Riemann Problem (GRP) of the diffusion equation. 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 [34,36,57–59], which leads to a sparse four-point block system for the scalar pressure. Once the new pressure field is known, the velocity vector field can subsequently be updated directly. A Picard iteration procedure that embraces the entire scheme in each time step is then used in order to achieve arbitrary high-order of accuracy also in time for the nonlinear convective term, without introducing a non linearity in the system for the pressure. The rest of the paper is organized as follows: in Section 2 the numerical method is described in detail, while in Section 4 a set of numerical test problems is solved in order to study the spatial and temporal accuracy of the presented approach. Some concluding remarks are given in Section 5. 2. Staggered space–time DG scheme for the 2D incompressible Navier–Stokes equations 2.1. Governing equations The two dimensional incompressible Navier–Stokes equations are given by

@v þ r  Fc þ rp ¼ mDv þ S; @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; u and v are the velocity components in the x and y direction, respectively; S ¼ Sðv ; x; y; tÞ is a (nonlinear) algebraic source term; Fc ¼ v  v is the flux tensor of the nonlinear convective terms, namely:

Fc ¼



uu

uv

v u vv

 :

Following the same idea of [56,60], the viscosity term is first written as mDv ¼ r  ðmrv Þ and then grouped with the nonlinear convective term. In this way the momentum Eq. (1) can be rewritten as

@v þ r  F þ rp ¼ S; @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. [56,60]. 2.2. Staggered unstructured grid Through this paper we use the same unstructured staggered grid in space as the one used in [36,37]. In the following, we briefly summarize the grid construction and the main notation for the spatial grid. After that, the primary and dual spatial elements are extended to the primary and dual space–time control volumes, respectively. 2.2.1. Unstructured staggered grid in space The spatial computational domain is covered with a set of N e non-overlapping triangular elements T i with i ¼ 1 . . . N e . By denoting with N d the total number of edges, the physical j-th edge will be called Cj . BðXÞ denotes the set of indices j corresponding to

237

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

boundary edges. The three edges of each triangle T i constitute the set Si defined by Si ¼ fj 2 ½1; N d  j Cj is an edge of T i g. For every j 2 ½1 . . . N d   BðXÞ there exist two triangles i1 and i2 that share Cj . We assign arbitrarily a left and a right triangle called ‘ðjÞ and rðjÞ, respectively. The standard positive direction is assumed to nj denote the unit normal vector defined be from left to right. Let ~ 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 that share the edge Cj is denoted by }ði; jÞ. For every j 2 ½1; N d   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 [36,57,59]. We denote by T i;j ¼ Rj \ T i the intersection element for every i and j 2 Si . Fig. 1 summarizes the used notation, the primal triangular mesh and the dual quadrilateral grid. According to [37], we will call the mesh of triangular elements fT i gi2½1;Ne  the main grid or primal grid and the quadrilateral grid fRj gj2½1;Nd  is termed the dual grid. The dual grid is covering X with non-overlapping quadrilaterals, so we define the equivalent quantities given for the main grid also for the dual one, briefly: N l is the total number of edges of Rj ; Cl indicates the physical l-th 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; ~ nl is the standard normal vector defined on l and assumed positive with respect to the standard orientation on l (defined, as for the main grid, 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 the point ‘ðjÞ and oriented in counter-clockwise direction. 2.2.2. Space–time extension In the time direction we cover the time interval ½0; T with a sequence of times 0 ¼ t0 < t1 < t2    < t N < t Nþ1 ¼ T. We denote the time step by Dt nþ1 ¼ tnþ1  tn and the corresponding time interval by T nþ1 ¼ ½t n ; t nþ1  for n ¼ 0 . . . N. In order to ease notation, sometimes we will use the abbreviation Dt ¼ Dtnþ1 . In this way the generic space–time element defined in the time interval ½t n ; t nþ1  is nþ1 nþ1 for the main grid, and Rst for given by T st i ¼ Ti  T j ¼ Rj  T the dual grid. Fig. 2 shows a graphical representation of the primary and dual space–time control volumes.

2.3. Space–time basis functions According to [36,37] we proceed as follows: we first construct the polynomial basis up to a generic polynomial degree p on some triangular and quadrilateral reference elements. In order to do this we take T std ¼ fðn; cÞ 2 R2;þ j c  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Þ 2

basis functions

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.

st Fig. 2. Example of space–time elements T st i (red) and Rj (green) with j 2 Si . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

interpolation polynomials passing through the Gauss–Legendre quadrature points for the unit interval. For every time interval ½tn ; tnþ1 , the map between the reference interval and the physical one is simply given by t ¼ t n þ sDtnþ1 for every s 2 ½0; 1. Using the tensor product we can finally construct the basis functions on the space–time elements T st and Rst such as i j ~ ~ /ðn; c; sÞ ¼ /ðn; cÞ  cðsÞ and wðn; c; sÞ ¼ wðn; cÞ  cðsÞ. The total st number of basis functions becomes N st / ¼ N /  N c and N w ¼ N w  N c .

2.4. Semi-implicit space–time DG scheme The discrete pressure ph is defined on the main grid, namely ph ðx; y; tÞjT st ¼ pi ðx; y; tÞ, while the discrete velocity vector field v h i

is defined on the dual grid, namely

vh ðx; y; tÞjR

st j

¼ v j ðx; y; tÞ.

f/k gk2½1;N/  on T std and N w ¼ ðp þ 1Þ2 basis functions on Rstd . The

The numerical solution of (2) and (3) is represented inside the space–time control volumes of the primal and the dual grid during

connection between reference and physical space is performed by the maps T i : T i ! T std for every i ¼ 1 . . . N e ; T j : Rj ! Rstd for

the current time interval T nþ1 by piecewise space–time polynomials as follows:

: Ti every j ¼ 1 . . . N d and its inverse, called T 1 i

T std and

: Rj Rstd , respectively. The maps from the physical coordiT 1 j nates to the reference one can be constructed following a classical sub-parametric or a complete iso-parametric approach. In the same way we construct the time basis functions on a reference interval ½0; 1 for polynomials of degree pc . In this case the resulting N c ¼ pc þ 1 basis functions fck gk2½1;Nc  are defined as the Lagrange

st

pi ðx; y; tÞ ¼

N/ X ~ ðiÞ ðx; y; tÞp ~ ðiÞ ðx; y; tÞp ^ nþ1 ^nþ1 / ¼: / ; i l;i l

ð4Þ

l¼1 st

vj ðx; y; tÞ ¼

Nw X ~ ðjÞ ðx; y; tÞv ~ ðjÞ ðx; y; tÞv ^ nþ1 ^ nþ1 w ¼: w ; j l;j l l¼1

ð5Þ

238

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

~ y; tÞ and wðx; ~ y; tÞ are genwhere the vectors of basis functions /ðx; ~ c; sÞ on T std  ½0; 1 and wðn; c; sÞ erated via the mappings from /ðn; on Rstd  ½0; 1, respectively. A weak formulation of the continuity Eq. (2) is obtained by mul~ and integrating over a control volume T st , for every tiplying it by / i st k ¼ 1 . . . N / . The resulting weak formulation for the discrete velocity v h reads

Z

T st i

~ ðiÞ / k

r  v h dxdydt ¼ 0:

ð6Þ

nþ1 nþ1 nij ¼ ~ ni jCst ; T st where ~ ; and Cst . We stress j ¼ Cj  T i;j ¼ T i;j  T j

again that the same weak form of the momentum Eq. (10) could have been also derived by forward and backward integration by parts of (7), following the ideas of Bassi and Rebay [68]. Using definitions (4) and (5), we conveniently rewrite the above equations as

X Z

  Z ~ ðjÞ @ v h þ r  Fh dxdydt þ ~ ðjÞ rp dxdydt w w h k k @t Rst Rst j j Z ~ ðjÞ Sh dxdydt; w ¼ k Rst j

~ ðiÞ v h  ~ / ni dsdt  k

Z T st i

~ ðiÞ  v h dxdydt ¼ 0; r/ k

ð8Þ

where ~ ni indicates the outward pointing unit normal vector. Since ph is allowed to jump at the element boundaries of the primary triangular grid and v h may jump at the boundaries of the dual quadrilateral elements, the integrals appearing in Eqs. (7) and (8) have to be split into elementwise contributions. Note, however, that thanks to the use of a staggered grid we do not need a Riemann solver here, since all the quantities are readily defined where needed for the flux computation. In particular, the velocity is continuous across the boundaries of the triangles on the main grid and the pressure is continuous across the boundaries of the dual quadrilateral grid. If we want to integrate the second term in (7) directly, without forward and backward integration by parts, the gradient of the discrete pressure ph appearing in the weak form of the momentum Eq. (7) can be easily interpreted in the sense of distributions, similar to a non-conservative product, since the discrete pressure is allowed 2 to jump at Cst j inside a dual velocity control volume. At this point we refer the reader to the extensive literature on high order accurate path-conservative finite volume and discontinuous Galerkin finite element schemes, where weak formulations of PDE with non-conservative products have already been derived and discussed in great detail, see e.g. [36,61–67]. We therefore obtain the following weak form of the continuity and the momentum equation:

X Z Cst j

j2Si

~ ðiÞ v j  ~ / nij dsdt  k



ð11Þ

Z

Cst j

Z

T st i;j

~ ðiÞ  v j dxdydt r/ k

¼ 0;

Rst j

T st rðjÞ;j

Rst j

~ ðjÞ S dxdydt; w k

Z Rst j

Cst j

~ ðjÞ Sdxdydt; w k

ð12Þ

Rj

# ðjÞ ~ wk v j dxdy

Z



t¼t nþ1

"Z Rj

# ðjÞ ~ wk v j dxdy t¼t n

ðjÞ

Rst j

~ @w k @t

v j dxdydt:

ð13Þ

In Eq. (13) we can recognize the fluxes between the current space– time element Rj  T nþ1 , the future space–time slab and the past space–time elements, as well as an internal space–time volume contribution that connects the layers inside the space–time element Rst j in an asymmetric way. Note how the asymmetry affects only the space–time volume contribution in (13). This is due to the nature of the time derivative operator, which has a natural positive direction given by the causality principle in time. Also note that the surface integral over the element at the past time includes the initial condition of the velocity in a weak sense. It can also be interpreted as using an upwind flux in time direction. By substituting Eq. (13) into (12) we obtain the following weak formulation of the momentum equation in space–time:

"Z Rj

# ~ ðjÞ w ~ ðjÞ dxdy w k l

þ þ ¼

Z

"Z



ð9Þ

  Z Z ~ ðjÞ @ v j þ r  Fj dxdydt þ ~ ðjÞ rp dxdydt w w ‘ðjÞ k k st @t Rst T j ‘ðjÞ;j Z Z   ~ ðjÞ rp dxdydt þ ~ ðjÞ p  p ~ þ w w rðjÞ rðjÞ ‘ðjÞ nj dsdt k k

~ ðjÞ / ~ ð‘ðjÞÞ~ ^nþ1 w nj dsdt p l;‘ðjÞ ¼ k l

~ ðjÞ @ v j dxdydt ¼ w k @t

!

Z

T st rðjÞ;j

where we have used the standard summation convention for the repeated index l. Integrating the first integral in (12) by parts in time and assuming the mesh to be fixed in time (i.e. the elements do not move), we obtain

and

2

¼ 0;

Z Z ~ ðjÞ @ v j dxdydt þ ~ ðjÞ r  Fdxdydt þ ~ ðjÞ r/ ~ ð‘ðjÞÞ dxdydt p ^nþ1 w w w l;‘ðjÞ k k k l st st @t Rst R T j j ‘ðjÞ;j Z Z ~ ðjÞ r/ ~ ðrðjÞÞ dxdydt p ~ ðjÞ / ~ ðrðjÞÞ~ ^nþ1 ^nþ1 w w þ nj dsdt p l;rðjÞ þ l;rðjÞ k l k l

ð7Þ

for every j ¼ 1 . . . Nd and k ¼ 1 . . . Nst w . Using integration by parts Eq. (6) becomes

¼

T st i;j

~ ðiÞ w ~ ðjÞ dxdydt v ^ nþ1 r/ l;j k l

Z

Z

@T st i

Cst j

j2Si

!

Z

and

~ and inteSimilarly, multiplication of the momentum Eq. (3) by w yields grating over a control volume Rst j

I

~ ðiÞ w ~ ðjÞ~ ^ nþ1  / k l nij dsdt v l;j

Z T st ‘ðjÞ;j

Z

"Z

Cst j

Rj

þ

! ~ ðjÞ ðjÞ @w k ~ ^ nþ1 wl dxdydt v l;j st @t R j t¼tnþ1 Z ~ ðjÞ r/ ~ ð‘ðjÞÞ dxdydt p ~ ðjÞ r/ ~ ðrðjÞÞ dxdydt p ^nþ1 ^nþ1 w w l;‘ðjÞ þ l;rðjÞ k l k l 

Z

~ ðjÞ / ~ ðrðjÞÞ~ ^nþ1 w nj dsdt p l;rðjÞ  k l #

~ ðjÞ w ~ ðjÞ dxdydt w k l

Z

Rst j

v^ nl;j  t¼t n

T st rðjÞ;j

Z Cst j

Z Rst j

~ ðjÞ / ~ ð‘ðjÞÞ~ ^nþ1 w nj dsdt p l;‘ðjÞ k l ~ ðjÞ r  F dxdydt w k

~ ðjÞ S dxdydt: w k

ð14Þ

For every i and j, Eqs. (11)–(14) are written in a compact matrix form as

Cst j

ð10Þ

Note that the term rph is not a true non-conservative product, since it can be written in a conservative divergence form as rph ¼ r  ðph IÞ, with the identity matrix I. Nevertheless, if ph contains discontinuities, the term rph can be properly interpreted as a special case within the more general framework of PDE with non-conservative products.

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

ð15Þ

j2Si

and

  ^ nþ1 ^ nþ1 ^ nj þ !j ðv ; rv Þ þ Rj p ^ nþ1 M þj  M j v  M j v j rðjÞ  Lj p‘ðjÞ ¼ Sj ; ð16Þ

239

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

respectively, where:

Z

M þj ¼

Rj

Z

M j ¼ M j ¼

!j ¼

Lj ¼

Sj ¼

Z

ð17Þ

~ ðjÞ ðx; y; tnþ1 Þw ~ ðjÞ ðx; y; t n Þdxdy; w k l

ð18Þ

ðjÞ

Rst j

Z

Rst j

Di;j ¼

Rj ¼

Rj

~ ðjÞ ðx; y; tnþ1 Þw ~ ðjÞ ðx; y; t nþ1 Þdxdy; w k l

Cst j

Cst j

Z Cst j

ð19Þ

~ ðjÞ r  F dxdydt; w k

Z

Z

~ @w k ~ ðjÞ w dxdydt; @t l

~ ðiÞ w ~ ðjÞ~ / k l nij dsdt 

~ ðjÞ / ~ ðrðjÞÞ~ w nj dsdt þ k l

~ ðjÞ / ~ ð‘ðjÞÞ~ w nj dsdt  k l

ð20Þ Z

~ ðiÞ w ~ ðjÞ dxdydt; r/ k l

T st i;j

Z T st rðjÞ;j

Z T st ‘ðjÞ;j

ð21Þ

~ ðjÞ r/ ~ ðrðjÞÞ dxdydt; w k l

ð22Þ

~ ðjÞ r/ ~ ð‘ðjÞÞ dxdydt; w k l

ð23Þ

Rst j

~ ðjÞ S dxdydt: w k

ð24Þ

Remark how introduces, for polynomial degrees pc > 0, an asymmetric contribution in time. The action of matrices L and R can be generalized by introducing a new matrix Qi;j , defined as

ri;j

Z

T st i;j

~ ðjÞ r/ ~ ðiÞ dxdydt  w k l

Z

X X ^ nþ1;kþ1 þ ^ nþ1;kþ1 Di;j M 1 Di;j M 1 j Qi;j pi j Q}ði;jÞ;j p}ði;jÞ X 1 2 ¼ Di;j Fc v nþ1;kþ ; j

j2Si

ð32Þ

Cst j

or, by introducing the boundary elements (see e.g. [37]),

"

X

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

j2Si \BðXÞ

X



X

ð25Þ

¼

j2Si BðXÞ

X

j2Si BðXÞ

ri;j is a sign function defined by

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

ð26Þ

  ^ nþ1 ^ nþ1 ^ nj þ !j ðv ; rv Þ þ QrðjÞ;j p ^ nþ1 M þj  M j v  M j v j rðjÞ þ Q‘ðjÞ;j p‘ðjÞ ¼ Sj ;

ð27Þ or, equivalently,

# ^ nþ1;kþ1 Di;j M 1 j Qi;j pi

^ nþ1;kþ1 Di;j M 1 j Q}ði;jÞ;j p}ði;jÞ

j2Si BðXÞ

~ ðjÞ / ~ ðiÞ ri;j~ w nj dsdt; k l

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

1

2 Di;j Fc v nþ1;kþ þ j

X

1

2 D@i;j Fc v nþ1;kþ ; j

ð33Þ

j2Si \BðXÞ

where D@i;j and Q@i;j are the natural extension of D and Q on triangular dual boundary elements, see e.g. [37]. Now the right hand side of Eq. (31) can be computed by using the velocity field at the old Picard iteration k and including the viscous effects using a fractional step type procedure. In this way, Eq. (31) represents a block ^ inþ1;kþ1 . Once the new presfour-point system for the new pressure p sure field is known, the velocity vector field at the new Picard iter^ nþ1;kþ1 can be readily updated from the momentum Eq. (30). ation v 2.5. Nonlinear convection–diffusion

  ^ nþ1 ^ nþ1 ^ nj þ !j ðv ; rv Þ þ Qi;j p ^ nþ1 M þj  M j v  M j v þ Q}ði;jÞ;j p j i }ði;jÞ ¼ Sj : ð28Þ In order to further ease notation, we will use the abbreviation  Mj ¼ Mþ j  M j henceforth and will write Eqs. (15) and (16) as follows:

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

ð29Þ

j2Si

^ nþ1 ^ nþ1 Mj v  M j Fc v j þ QrðjÞ;j p^ nþ1 j rðjÞ þ Q‘ðjÞ;j p‘ðjÞ ¼ 0;

ð30Þ

where Fc v j is an appropriate discretization of the nonlinear convective, viscous and source terms. The details for the computation of v j will be presented later. Formal substitution of the discrete Fc momentum Eq. (30) into the discrete continuity Eq. (29), see also [35–37,39], yields

X X X ^ nþ1 þ ^ nþ1 Di;j M 1 Di;j M 1 Di;j Fc v j: j Qi;j pi j Q}ði;jÞ;j p}ði;jÞ ¼ j2Si

implicitly, then system (31) becomes nonlinear and it would be very cumbersome to solve it. In order to overcome this problem we introduce a simple Picard iteration to introduce the information of the new pressure into the viscous and convective terms, but without introducing a nonlinearity in the final system to be solved. This approach is inspired by the local space–time Galerkin predictor method proposed for the high order time discretization of PN P M schemes in [60,69]. Hence, for k ¼ 1; Npic , we rewrite system (31) as

j2Si

M j

where

convective-viscous contribution in the time interval T nþ1 is based on the old information T n and does not see the effects of the new pressure in the time interval T nþ1 . Furthermore, if we take Fc vj

j2Si

Z

Qi;j ¼

at time tn and hence Eq. (31) would represent a four-point block ^ nþ1 system for the new pressure p , as proposed in [37]. i Unfortunately, in problems when the convective-viscous effects cannot be neglected, this will produce only a low order accurate method in time. The problem in this case is that the

j2Si

ð31Þ

j2Si

We have now to choose a time discretization for the nonlinear convective-viscous term. The simplest choice would be to take P Fc v j explicitly, so in this case Di;j Fc v n becomes a known term j2Si

j

To close the problem it remains to specify how to construct the 1 2 v nþ1;kþ . Following the nonlinear convection–diffusion operator Fc j ideas of [37], a space–time DG scheme for the convection–diffusion terms on the dual mesh is given by

Z

Z ~ k @ v h dxdydt þ ~ k Gh  ~ w w n dsdt @t Rst @Rst j j Z Z ~ k  Fðv h ; rv h Þdxdydt ¼ ~ ðjÞ S dxdydt; w rw  k Rst j

Rst j

ð34Þ

and the numerical flux for both, the convective and the viscous contribution, is given such as in [55,56,60], and reads

n¼ Gh  ~

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

ð35Þ

with

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

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

ð36Þ

which contains the maximum eigenvalue of the Jacobian matrix of the purely convective transport operator Fc in normal direction, see [35], and the stabilization term for the viscous flux, see [56,60].

240

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

Furthermore, the v h and rv h denote the velocity vectors and their gradients, extrapolated to the boundary of Rj from within the eleþ



ment 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. We discretize the velocity v h explicitly but its gradient has to be taken implicitly, in order to avoid additional restrictions on the maximum time step given by the viscous terms. In viscosity dominated problems, this allows us to use both, high viscosity and large time steps. After integration of the first term of (34) by parts in time the resulting fully discrete formulation of (34) becomes

v^ nþ1;kþ j

1 2

nþ1;kþ12

 n 1 nþ1;k ^ ¼ M 1 ; rv h j M j v j  M j !j ðv h

Þ þ M 1 j Sj ;

ð37Þ

where

!j ðv h ; rv h Þ ¼ ¼

Z Rst j

~ ðjÞ r  Fðv h ; rv h Þdxdydt w k

Z

@Rst j

~ k Gh  ~ w n dsdt 

Z Rst j

~ k  Fðv h ; rv h Þdxdydt: ð38Þ rw

Due to the explicit treatment of the nonlinear convective terms, 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

ð39Þ

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 [70,71]. Implicit discretization of the viscous contribution rv in (34) involves two five-point block systems (one for each velocity component) that can be efficiently solved using a matrix-free GMRES algorithm [72]. The solution of this system is not necessary in problems where the viscous term is small enough to be integrated explicitly in time. In that case, i.e. for explicit discretizations of the viscous terms, one has to include the additional explicit time step restriction for parabolic PDE in Eq. (39). 1 1 nþ1;kþ12 Once v has been computed, we set Fc v nþ1;kþ2 :¼ v^ nþ1;kþ2 . As j

j

^ nþ1;0 we can take the old velocity initial guess v j lation of

j

v nh , or the extrapo-

v nh into the interval T nþ1 .

2.6. Pressure correction formulation and final algorithm The preliminary algorithm described above, as formulated by Eqs. (37), (33) and (30) still contains an important drawback: indeed, Eq. (37) does not depend on the pressure of the previous Picard iteration and hence the algorithm does not see the effect of the pressure in the time interval T nþ1 . In order to overcome the problem we introduce the contribution of the pressure from the previous Picard iteration directly into Eq. (37). Then, we update ^ inþ1;kþ1 . With this modification, the velocity with the new pressure p Eqs. (37), (33) and (30) and hence the final algorithm become: 1 2

v^ nþ1;kþ j

nþ1;kþ12

 n 1 nþ1;k ^ ¼ M 1 ; rv h j M j v j  M j ! j ðv h

^ nþ1;k þ M 1  M 1 j Q‘ðjÞ;j p‘ðjÞ j Sj ;

^ nþ1;k Þ  M 1 j QrðjÞ;j prðjÞ ð40Þ

  X   X nþ1;k ^ nþ1;kþ1  p ^ inþ1;k þ ^ nþ1;kþ1  p ^ }ði;jÞ Di;j M 1 Di;j M 1 j Qi;j pi j Q}ði;jÞ;j p}ði;jÞ j2Si

¼

X j2Si

j2Si nþ1;kþ12

Di;j Fc vj

;

v^ nþ1;kþ1 ¼ Fc v nþ1;kþ j j

   i nþ1;kþ1 ^ nþ1;kþ1 ^ nþ1;k ^ ‘ðjÞ ^ nþ1;k þ Q‘ðjÞ;j p :  M 1 QrðjÞ;j p p p j rðjÞ rðjÞ ‘ðjÞ

ð42Þ

Note that Eqs. (41) and (42) are written in terms of the pressure cor  ^ inþ1;kþ1 ¼ p ^ nþ1;kþ1 ^ nþ1;k . Since both space and time are p rection Dp i i involved in this algorithm and since space and time discretization are intrinsically coupled, we call this method a coupled space–time pressure correction algorithm. The resulting linear system for the pressure correction is very sparse thanks to the use of the staggered grid, including only four non-zero blocks per element. Note that the same algorithm on a collocated grid would lead to a pressure system with 10 non-zero blocks per element, since substituting the discrete velocity into the continuity equation on a collocated grid involves also neighbors of neighbors. The significantly improved sparsity pattern of the linear system is indeed a key advantage of the algorithm presented in this paper. In all subsequent numerical examples, we were able to solve the system (41) using a simple matrix-free version of the GMRES algorithm [72] without using any preconditioner, which is a unique property in the context of high order implicit DG schemes. As initial guess for the pressure one can take pnþ1;0 ¼ 0, but one h could also choose the extrapolation of pnh into T nþ1 . One time step of the final algorithm can be summarized as follows: 1. Initialize v nþ1;0 and pnþ1;0 using the known information from T n ; h h 2. Picard iteration over k ¼ 0 . . . N pic : nþ1;kþ1

2 using (40), i.e. convective terms are dis(a) compute v h cretized explicitly and viscous terms implicitly; then set 1 1 v nþ1;kþ2 :¼ v^ nþ1;kþ2 , Fc

j

j

^ nþ1;kþ1 by solving the discrete pressure Poisson (b) compute p Eq. (41), ^ jnþ1;kþ1 explicitly from (42); (c) update v ^ nþ1 ¼ p ^ nþ1;kþ1 . ^ nþ1 ^ jnþ1;kþ1 and p ¼v 3. set v j For the spatial computational domain we can apply the remark given in [37] and so either use a subparametric or a complete isoparametric representation. The second approach requires to store more information about each element, but it also allows to generalize the shape of the elements. This property is crucial when we try to discretize complex curved domains with a very coarse grid. In any case, this generalization does not affect the computational time during run-time, since it interests only the construction of the geometry-dependent matrices in the preprocessing stage of the algorithm. 2.7. Splitting of the space–time matrices into a spatial and temporal part Even if the shape of the main matrices is similar compared to the ones introduced in [37], the number of degree of freedom and the integral values are, in general, different. Due to the tensor product construction of the space–time basis functions, we can split the main integrals (17)–(21) and (25) into a spatial and a temporal part. Briefly, the space–time matrices are generated from the spatial matrices of [37], componentwise, as:

M þj ðk; lÞ ¼ c‘2 ðkÞ ðt nþ1 Þc‘2 ðkÞ ðtnþ1 ÞM sj ð‘1 ðkÞ; ‘1 ðlÞÞ;

ð43Þ

M j ðk; lÞ ¼ c‘2 ðkÞ ðt nþ1 Þc‘2 ðkÞ ðtn ÞM sj ð‘1 ðkÞ; ‘1 ðlÞÞ;

ð44Þ

M j ðk; lÞ ¼ M sj ð‘1 ðkÞ; ‘1 ðlÞÞDt ð‘2 ðkÞ; ‘2 ðlÞÞ;

ð45Þ

Di;j ðk; lÞ ¼ Dt Dsi;j ð‘1 ðkÞ; ‘1 ðlÞÞM t ð‘2 ðkÞ; ‘2 ðlÞÞ; Qi;j ðk; lÞ ¼ Dtnþ1 Qsi;j ð‘1 ðkÞ; ‘1 ðlÞÞM t ð‘2 ðkÞ; ‘2 ðlÞÞ;

ð47Þ

nþ1

ð41Þ

1 2

h

ð46Þ

241

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

where the apex s means that the matrix is the one constructed in [37]; Dt and M t are two time matrices defined as

  Z ~ ~l ¼ Dt k;   ~ ~l ¼ M t k;

1

0

Z

dck~ ðnÞ c~l ðnÞdn; dn

ð48Þ

ck~ ðnÞc~l ðnÞdn;

ð49Þ

1

0

and ‘1 ; ‘2 are two appropriate sorting functions defined according to the number of space and time basis functions. Remark how the action of the matrix Dt defined in (48) is symmetric only if pc ¼ 0. 3. Stability analysis of the space–time DG method In this section we present a proof regarding the stability in the space–time DG framework. In the particular case of pc ¼ 0, the stability follows directly from the stability proof given in [37]. In the general case of pc > 0 we can derive the following result:

guess is taken according to (52), i.e. p0h ¼ 0. We can further suppose that S ¼ 0 and m ¼ 0 so that Fh ðv h ; rv h Þ ¼ Fh ðv h Þ since the source term could increase naturally the total kinetic energy of the system and we want a stability result that is independent of the source term or the physical viscosity. From (51) and (52) and taking into account the upwinding in time direction, we get

X Z Cst j

j2Si

~ ðiÞ v h  ~ / nij ds dt  k

!

Z T st i;j

~ ðiÞ  v h dx dt r/ k

¼ 0;

ð50Þ

#

Rj

~ ðjÞ v h dx w k

þ þ 

 t¼t nþ1

Z

Cst j

Z

Rst j

Rj

# ~ ðjÞ v h dx w k

~ ðjÞ rp dx dt þ w ‘ðjÞ k

T st ‘ðjÞ;j

Z

"Z

t¼t n

Z T st rðjÞ;j

Z Rst j

ðjÞ

Rst j

~ @w k @t

vh dx dt

~ ðjÞ rp dx dt w rðjÞ k

Z   ~ ðjÞ p  p ~ w n ds dt þ j rðjÞ ‘ðjÞ k

~ ðjÞ  Fh dx dt ¼ rw k



Z

@Rst j

~ ðjÞ Gh  n ds dt w k

~ ðjÞ Sh dx dt: w k

ð51Þ

where we have used the abbreviation dx ¼ dxdy. In Eq. (50) we proceed as in [37] and take ph as test function. Passing to the sum over each element T st j and exchanging the two summations we obtain

X Z



Cst j

j2Si

ph;i v h  ~ nij ds dt 

Ne X Z X

Cst j

i¼1 j2Si

T st ‘ðjÞ;j

þ

Z Cst j

ph;i v h  ~ nij ds dt 

Z T st i;j

rph;i  v h dx dt

rph;‘ðjÞ  v h dx dt þ

Z T st rðjÞ;j

@Rst j

@vh 1  v h dx dt ¼ 2 @t Rst j "Z # 1  v þh  vþh dx 2 Rj

v v þ h 

 h dx t¼tn

v h  ðGh  nÞdsdt 

!

Z Rst j

rv h : Fh dxdt ¼ 0; ð53Þ

Z Rst j

@ 2 1 v dx dt ¼ @t h 2

"Z Rj

#

v h  v h dx t¼tnþ1

ð54Þ

:

t¼t n

In this way, Eq. (53) becomes "Z

"Z # 1   v v  v v  v  v dx 2 Rj h h Rj Rj j¼1 t¼t nþ1 t¼t n t¼t nþ1 "Z # ! Z Z 1 v þh  v þh dx þ v rv h : Fh dxdt ¼ 0: ð55Þ þ h  ðGh  nÞdsdt  2 Rj @Rst Rst n j j

Nd X

#

 h

"Z

 h dx

#

þ h

 h dx

t¼t

"Z # "Z # "Z # Nd X 1 1 1 v h vh dx  v h v h dx þ v h v h dx 2 Rj 2 Rj 2 Rj j¼1 t¼tnþ1 t¼t n t¼tn "Z # "Z # 1 vþ vþ dx  v þh vh dx þ 2 Rj h h Rj t¼t n t¼tn ! Z Z ð56Þ þ v h  ðGh  nÞ ds dt  rv h : Fh dx dt ¼ 0: @Rst j

Rst j

We can recognize that

þ

"Z

1 2

Rj

1 ¼ 2

#

v h v h dx

þ t¼t n

"Z

 Rj

vþh  v h

2

"Z # "Z # 1 v þh v þh dx  v þh v h dx 2 Rj Rj t¼tn t¼tn #

dx

P 0;

ð57Þ

t¼t n

and hence

+

!

Rst j

Z

@vh  v h dxdt þ @t

# Rj

t¼t nþ1

Z

¼ 0;

¼ 0; m

" Nd Z X j¼1

T st i;j

rph;i  v h dx dt

Z



"Z

þ where v  h jt¼t and v h jt¼t are the left and the right value of v h at the nt , respectively. We interface according to the time normal vector ~ can now write the third integral in (53) as

!

Z

v v

 h dx

Following [73] we rearrange Eq. (55) by adding and subtracting the R  quantity 12 v  h v h dx t¼t n and thus obtain

and

"Z

#  h

Rj

j¼1

Theorem 1. The proposed staggered space–time DG method (8) and (14) is L2 stable if the discretization of the nonlinear convective and viscous terms is L2 stable. Proof 1. By construction, the numerical solution satisfies the weak formulation of the continuity and momentum equation, namely:

"Z

Nd X

ð52Þ

rph;rðjÞ  v h dx dt

#   ph;rðjÞ  ph;‘ðjÞ v h  ~ nj ds dt ¼ 0:

We now take v h as test function in (51) and sum up both components of the momentum equation. Furthermore, we sum up over all dual elements. Then, one obtains thanks to (52) that all the pressure contributions add up to zero, as in [37]. Note that this property is also satisfied for the pressure at each Picard iteration, if the initial

"Z # "Z # Nd X 1 1     v v dx  v v dx 2 Rj h h 2 Rj h h j¼1 t¼tnþ1 t¼t n ! Z Z v h  ðGh  nÞ ds dt  rv h : Fh dx dt 6 0: þ @Rst j

Rst j

ð58Þ

One can observe how in (58) the difference of the kinetic energy is computed at the end of two time slices. As for the method proposed in [37], the stability of the scheme is therefore determined by the stability of the discretization used for the nonlinear convective term, which is assured according to [74,75] when appropriate monotone entropy fluxes are used. One thus obtains the result:

1 2

Z X



vh vh dx

 t¼t nþ1

1 2

Z X

v h v h dx

 6 0:  t¼tn

ð59Þ

242

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

Table 1 Numerical convergence results for the manufactured solution test problem with polynomial degrees p ¼ 1 to p ¼ 4 in space and time. Ne

ðpÞ

40 160 640 2560

1.217E01 2.678E02 6.050E03 1.758E03

40 160 640

7.703E04 3.864E05 2.425E06

ðv Þ

rðpÞ

rðv Þ

9.572E02 2.362E02 5.527E03 1.497E03

– 2.2 2.1 1.8

– 2.0 2.1 1.9

1.425E03 4.999E05 1.974E06

– 4.3 4.0

– 4.8 4.7

ðpÞ

p ¼ pc ¼ 1

ðv Þ

rðpÞ

rðv Þ

1.052E02 1.065E03 9.103E05 7.820E06

– 3.3 3.1 3.0

– 3.3 3.5 3.5

7.135E05 1.418E06 2.945E08

– 5.5 5.2

– 5.7 5.6

p ¼ pc ¼ 2 8.740E03 8.833E04 1.050E04 1.347E05

p ¼ pc ¼ 3

p ¼ pc ¼ 4 5.315E05 1.143E06 3.102E08

Finally, note that even for Gh ¼ Fh 0 the high order space–time DG scheme in general dissipates the total kinetic energy, since the solution is allowed to jump at the two sides of a time slice and hence, in þ general, v  h jtn – v h jt n .

u0 ¼ v 0 ¼ 1; p0 ¼ 1; x ¼ 2p; k ¼ 10=2p; t end ¼ 0:5; Dt according to condition (39); and m ¼ 0:01. The temporal accuracy is chosen equal to the spatial one, the total number of Picard iterations is taken as Npic ¼ p þ 1 and pnþ1;0 0 for the present test. The computational

4. Numerical test problems

domain is X ¼ ½0:5; 0:52 ; the exact velocity field and pressure are taken as initial conditions and the exact pressure is also specified on @ X as boundary condition. The L2 error between the analytical and the numerical solution is computed as

In this section we study the accuracy of our new numerical method by solving some classical numerical benchmark problems, such as the lid-driven cavity flow, the unsteady oscillatory flow in a pipe or the unsteady flow past a circular cylinder. In particular, we perform quantitative comparisons between the numerical solution and available exact analytical solutions wherever possible. 4.1. Convergence test using a manufactured solution In order to study the accuracy of the proposed space–time DG method, we need an exact unsteady solution of (1)–(3). For that purpose, we propose a so-called manufactured solution in this section, which also makes use of a linear source term of the type Sðx; y; tÞ. The exact analytical solution for the velocity and the pressure is constructed so that

v an ¼ v 0 sin ½kðx  yÞ  xt;

pan ¼ p0 sin ½kðx  yÞ  xt;

ð60Þ

with the amplitudes v 0 ¼ ðu0 ; v 0 Þ and p0 . Using the manufactured solution ðv an ; pan Þ we can compute all terms in (1) exactly and hence obtain a source term Sðx; y; tÞ that balances the momentum equation. Remark that the velocity field must be divergence-free (r  v ¼ 0), hence u0 ¼ v 0 . In the present test case, we take

ðpÞ ¼

sZ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðph  pan Þ2 dxdy;

ðvÞ ¼

X

sZ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðv h  v an Þ2 dxdy; X

where the subscript h refers to the numerical solution obtained at the final time t ¼ t end . The resulting rate of convergence is shown in Table 1. We observe that the optimal order of convergence is obtained for even order schemes up to p ¼ 4 for the present unsteady test. 4.2. The Womersley problem Here we consider an unsteady, viscosity-dominated test problem for which the incompressible Navier–Stokes equations have a nontrivial exact solution, namely the fluid flow inside a rigid planar pipe that is driven by a sinusoidal pressure gradient of the type

! e P pout ðtÞ  pin ðtÞ ixt ¼R e : L q

ð61Þ

e is the amplitude of the presIn this test L denotes the tube length; P sure oscillation; q is the density of the fluid; x is the frequency of

Table 2 Numerical convergence results for the planar Womersley problem. p

pc

Ne

Nt

ðpÞ

ðv Þ

rðpÞ

rðv Þ

1 1 1 1

1 1 1 1

46 184 736 2944

6 12 24 48

5.7880182E02 1.7635947E02 4.6206559E03 1.1683966E03

1.8848423E03 5.5901107E04 1.4587701E04 3.7404869E05

– 1.7 1.9 2.0

– 1.8 1.9 2.0

2 2 2 2

2 2 2 2

46 184 736 2944

6 12 24 48

7.0716231E03 4.8160864E04 3.0677533E05 1.9295385E06

2.6412698E04 3.8846170E05 7.2036760E06 1.6070616E06

– 3.9 4.0 4.0

– 2.8 2.4 2.2

3 3 3 3

3 3 3 3

46 184 736 2944

6 12 24 48

9.8372146E04 7.7144497E05 5.0814347E06 3.2173776E07

1.2793693E05 7.8462176E07 4.8795894E08 3.0326872E09

– 3.7 3.9 4.0

– 4.0 4.0 4.0

4 4 4 4

4 4 4 4

46 184 736 2944

6 12 24 48

7.3692980E05 1.2539784E06 2.1930727E08 1.0258845E09

5.1193160E07 2.1649081E08 1.1576584E09 7.0131498E11

– 5.9 5.8 4.4

– 4.6 4.2 4.0

243

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

0.25

1.5

DG (p=1) DG (p=2) DG (p=3) Reference

0.2

DG (p=1) DG (p=2) DG (p=3) Reference

1

0.15 0.1

0.5

0

p

u

0.05 0

-0.05 -0.5

-0.1 -0.15

-1

-0.2 -0.25

0

0.25

0.5

0.75

1

1.25

-1.5

1.5

0

0.25

0.5

0.75

t

1

1.25

1.5

t

Fig. 3. Time series for the axial velocity u and the pressure p computed at ðx; yÞ ¼ ð0:5; 0Þ for the coarsest grid N e ¼ 46 and N t ¼ 6.

0.2

DG (p3p3) x=-0.2 DG (p3p3) x=0.3 Exact Solution

6

0.15

5 0.1

4

0

y

y

0.05

3

-0.05

2

-0.1 -0.15

1 -0.2 -0.2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

u

0

Fig. 4. Radial velocity profiles for x ¼ 0:2 and x ¼ 0:3 at times, from left to right, t ¼ ½0:75; 0:5; 0:875; 1:0; 0:125. Comparison between exact and numerical solution.

0

1

2

3

4

5

6

x Fig. 5. Velocity field of the Taylor–Green vortex on the coarse grid N e ¼ 40 with p ¼ 4. The edge-based dual grid is shown.

the oscillation; pin and pout indicate the inlet and the outlet pressure, respectively; R is the real part operator. By imposing Eq. (61) at the tube ends, the exact analytical solution for the three dimensional, axially symmetric case was found by Womersley in [76]. It can be derived also for the two dimensional planar case. The resulting axial velocity is uniform in the x-direction and is given by

" uðx; y; tÞ ¼ R i

where k ¼ bottom.

e P

q



x 1

# cos½kðyc  1Þ ; cosðkÞ

ð62Þ

pffiffiffiffiffiffiffiffiffiffiffi pffiffiffi b ia2 ; a ¼ R xm; yc ¼ yy ; and yb is the y value of the R

For the present test X ¼ ½0:5; 0:5  ½0:2; 0:2; and

eP q

¼ 1. We

take a set of successively refined grids in order to show the convergence behavior to the exact solution with respect to the order p in space and pc in time. According to [76] the nonlinear convection effect is neglected for the present test. Thus, the stability of our scheme is not restricted by the CFL condition on the fluid velocity. Since we use a high viscosity coefficient in this test, the implicit

treatment of the viscous term is necessary to allow large time steps. In particular we choose m ¼ 5  102 and t end ¼ 1:5. On the coarsest grid we use Dt ¼ tend =6, then the time step is reduced proportional to the spatial grid size. No-slip boundary conditions are imposed on the top and the bottom boundary, while the pressure (61) is imposed at the inlet and the outlet boundary on the left and on the right, respectively. The number of Picard iterations is given by N p ¼ p þ 1 for all simulations. The resulting convergence results, using the L2 -norm as in the previous example, are shown in Table 2. Observe how a non-optimal order of convergence p is achieved for the velocity for odd order schemes (even polynomial degree p), while the optimal convergence rate p þ 1 is achieved for the pressure for all polynomial degrees. Note that when using the semi-implicit staggered DG method introduced in [36] only a second order of convergence could be achieved for this unsteady test problem, while full high order convergence in space and time is obtained with the new scheme presented in this paper. In Fig. 3 we show the time series

244

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

Table 3 Numerical convergence results for the velocity vector field of the Taylor–Green vortex. Ne

40 160 640 2560

p ¼ pc ¼ 1

p ¼ pc ¼ 2

p ¼ pc ¼ 3

p ¼ pc ¼ 4

ðv Þ

r

ðv Þ

r

ðv Þ

r

ðv Þ

r

3.088E01 8.868E02 2.267E02 5.476E03

– 1.8 2.0 2.0

5.588E02 5.765E03 7.052E04 8.452E05

– 3.3 3.0 3.1

5.895E03 4.730E04 2.387E05 1.312E06

– 3.6 4.3 4.2

1.669E03 3.109E05 6.233E07 1.297E08

– 5.7 5.6 5.6

Fig. 6. Vorticity pattern for the double shear layer test at times, from top left to bottom right, t ¼ 0:4; t ¼ 0:8; t ¼ 1:2; t ¼ 1:8.

of the axial velocity and the pressure in a given point for the coarsest grid configuration ðN e ; N t Þ ¼ ð46; 6Þ. While piecewise linear space–time polynomials are not able to reproduce the sinusoidal signal well with only six time steps, the piecewise quadratic and higher order approximations in space and time yield an almost perfect match with the exact solution even on this extremely coarse space–time grid. In Fig. 4 we compare the resulting numerical velocity profiles uðyÞ against the exact solution at several times for the case ðp; pc Þ ¼ ð3; 3Þ and N e ¼ 736. Two different locations, x ¼ 0:2 and x ¼ 0:3, are plotted in order to show that the profile is constant in the x-direction. One observes that there is no visible difference between numerical and exact solution in Fig. 4.

4.3. Taylor–Green vortex Another widely used testcase for the verification of numerical methods for the incompressible Navier–Stokes equations is the Taylor–Green vortex problem. The analytical unsteady solution is given by

uðx; y; tÞ ¼ sinðxÞ cosðyÞe2mt ;

v ðx; y; tÞ ¼  cosðxÞ sinðyÞe2mt ;

ð63Þ

ð64Þ 1 ð65Þ pðx; y; tÞ ¼ ðcosð2xÞ þ cosð2yÞÞe4mt : 4 2 The computational domain is X ¼ ½0; 2p and is extended using

periodic boundary conditions on all the boundaries.

245

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249 v Ref (Ghia) u Ref (Ghia) u (DG p3p3) v (DG p3p3)

Re 100 1.2

0.4 1 0.8

0.2

0

y

Velocities

0.6 0.4 0.2 -0.2 0 -0.2 -0.4

-0.4

-0.4

-0.2

0

0.2

-0.4

0.4

-0.2

0

0.2

0.4

0.2

0.4

x

xi v Ref (Ghia) u Ref (Ghia) u (DG p3p3) v (DG p3p3)

Re 400 1.2

0.4 1 0.8 0.2

0.4 0

y

Velocities

0.6

0.2 0

-0.2

-0.2 -0.4 -0.6

-0.4

-0.4

-0.2

0

0.2

-0.4

0.4

-0.2

0

x

xi

Fig. 7. Velocity profiles (left) and streamlines (right) at Reynolds numbers Re ¼ 100 and Re ¼ 400 for the lid-driven cavity problem.

As implied by Eqs. (63)–(65), the resulting velocity field initially appears as depicted in Fig. 5 and then starts to lose energy according to the friction effects. For the present test we consider several grid refinements; tend ¼ 0:1; m ¼ 0:1; and Dt is chosen according to the CFL time restriction for the nonlinear convective terms. The numerical convergence results are shown in Table 3. We find that the optimal convergence rates are achieved for this important nontrivial test problem with periodic boundary conditions. 4.4. Double shear layer The numerical scheme is applied here to a test case studied in [77], which contains a high initial velocity gradient. We take

X ¼ ½1; 12 and, as initial condition, we consider a perturbed double shear layer profile:

u0 ¼

~ ðyn  0:25Þ; if yn 6 0:5; tanh ½q ~ ð0:75  yn Þ; if yn > 0:5; tanh ½q

ð66Þ

v 0 ¼ d sinð2pxn Þ;

ð67Þ

p0 ¼ 1;

ð68Þ

where yn ¼ yþ1 and xn ¼ xþ1 are the normalized vertical and horizontal 2 2 ~ is a parameter that determines the slope coordinates, respectively; q of the shear layer; and d is the amplitude of the initial perturbation. ~ ¼ 30; m ¼ 2  104 ; p ¼ 4 and For the present test we set d ¼ 0:05; q pc ¼ 3. The time step is chosen according to the CFL condition for the nonlinear convective terms and four Picard iterations have been used in this simulation. The domain X is covered with a total number of only Ne ¼ 640 triangles and periodic boundary conditions are imposed everywhere. The resulting vorticity pattern is reported at several times in Fig. 6. The two thin shear layers evolve into several

246

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

vortices, as observed in [77], and overall the small flow structures seem to be relatively well resolved also at the final time t ¼ 1:8, even if a very coarse grid has been used in space and time.

4.5. Lid-driven cavity flow We consider here another classical benchmark problem for the incompressible Navier–Stokes equations, namely the lid-driven cavity problem [78]. This test case is solved numerically with the new staggered space–time DG scheme on very coarse grids using polynomial degrees of p ¼ 3 and pc ¼ 3 in space and time, respectively. 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. at y ¼ 0:5) and impose so-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 only N e ¼ 116 triangles for Re ¼ 100; 400; 1000 and N e ¼ 512 triangles for Re ¼ 3200.

Re 1000

For the present test Dt is taken according to the CFL condition (39) and tend ¼ 150. According to [78], primary and corner vortices appear from Re ¼ 100 to Re ¼ 3200, a comparison of the velocities against the data presented by Ghia et al. in [78], as well as the streamline plots are shown in Figs. 7 and 8. A very good agreement is obtained in all cases, even if a very coarse grid has been used in space and time. 4.6. Flow over a circular cylinder In this section we consider the flow over a circular cylinder. In this case, the use of an isoparametric finite element approach is mandatory to represent the curved geometry of the cylinder wall, see [36,68]. We consider here the viscous case in order to show the formation of the von Karman vortex street. We take a sufficiently pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi large domain X ¼ ½20; 80  ½20; 20  f x2 þ y2 6 1g and we cover it with only N e ¼ 1702 triangles. Note that the chosen grid is extremely coarse compared to the dimension of the domain X.

v Ref (Ghia) u Ref (Ghia) u (DG p3p3) v (DG p3p3)

1.2

0.4

1 0.8

0.2

y

0.4 0.2

0

0 -0.2

-0.2

-0.4 -0.6 -0.8

-0.4 -0.4

-0.2

0

0.2

0.4

-0.4

-0.2

0

xi

0.2

0.4

x

Re 3200

v Ref (Ghia) u Ref (Ghia) u (DG p3p3) v (DG p3p3)

0.4

0.2

y

Velocities

0.6

0

-0.2

-0.4

.4

-0.2

0

xi

0.2

0.4

-0.4

-0.2

0

0.2

0.4

x

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

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

247

Fig. 9. Streamlines along the circular cylinder at times, from top left to bottom right, t ¼ 25; 50; 100 and t ¼ 200.

The characteristic average size of the mesh is h ¼ 1:295 and the smallest element size is about hmin ¼ 0:347. As initial condition  ; 0Þ, where u  is the inlet velocity, taking we set v ðx; y; 0Þ ¼ ðu  ¼ 0:5 in our case. For the present test we use Dt according to u  ; 0Þ is prescribed at the left (39); p ¼ 3; pc ¼ 2. The velocity ðu boundary while passage boundary conditions are imposed on the other external edges of the domain. Finally a no-slip wall boundary condition is imposed on the cylinder surface. A plot of the streamlines is reported in Fig. 9 at several output times. The resulting profiles for the vorticity and the horizontal velocity u are plotted in Fig. 10, as well as the dual grid elements for Re ¼ 100. As shown in Fig. 9, two vortices are initially generated at the circular cylinder and then, several vortices leave the cylinder and generate the Von Karman street as we can see in Fig. 10. The resulting Strouhal number for the present test is St ¼ fdu ¼ 0:1647 that is in good agreement with St ¼ 0:1649 obtained by Qu et al. in [79].

5. Conclusions A novel high order accurate staggered semi-implicit space–time discontinuous Galerkin scheme has been proposed for the solution of the two-dimensional incompressible Navier–Stokes equations

on unstructured curved triangular meshes. The use of a staggered grid makes our scheme different from the space–time DG schemes proposed in [16,17]. The high order in space and time was verified up to p ¼ 4 against available exact solutions for several test cases that include a manufactured solution using source terms, the viscosity-dominated Womersley problem and the well-known Taylor–Green vortex problem with periodic boundary conditions. The numerical results agree very well with the reference data for all test cases under consideration. In the special case pc ¼ 0 the numerical method proposed in this paper reduces exactly to the semi-implicit staggered DG scheme forwarded in [37], so it can be seen as its natural extension to high order of accuracy in time. Furthermore, the use of matrices that depend only on the geometry and on the polynomial degree and that hence can be precomputed before runtime, as well as a very good sparsity pattern involved in the solution of the main system for the pressure, leads to a computationally efficient scheme. The use of a staggered mesh leads to only four non-zero blocks per elements, while the same algorithm on a collocated mesh would result in 10 non-zero blocks per element, since also neighbors of neighbors would be involved. Actually, we have solved all our test problems with a matrix-free implementation of the GMRES method [72], without the use of any preconditioner, which is a very interesting feature of our

248

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249

10

y

5

0

-5

-10

-5

0

5

10

15

20

25

30

35

40

45

25

30

35

40

45

x 10

y

5

0

-5

-10

-5

0

5

10

15

20

x Fig. 10. Laminar viscous flow past a circular cylinder. Profile for the vorticity and horizontal velocity u at time t ¼ 300.

staggered space–time DG scheme compared to other implicit and semi-implicit DG schemes that exist in the literature. Future research will concern the extension of the scheme to the fully three-dimensional case on unstructured tetrahedral and hexahedral meshes and its application to turbulent flows. Acknowledgment The research presented in this paper was partially 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. References [1] Harlow F, Welch J. Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface. Phys Fluids 1965;8:2182–9. [2] Patankar V, Spalding B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int J Heat Mass Transfer 1972;15:1787–806. [3] Patankar V. Numerical heat transfer and fluid flow. Hemisphere Publishing Corporation; 1980. [4] van Kan J. A second-order accurate pressure correction method for viscous incompressible flow. SIAM J Sci Stat Comput 1986;7:870–91. [5] Taylor C, Hood P. A numerical solution of the Navier–Stokes equations using the finite element technique. Comput Fluids 1973;1:73–100. [6] Brooks A, Hughes T. Stream-line upwind/Petrov Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equation. Comput Methods Appl Mech Eng 1982;32:199–259. [7] Hughes T, Mallet M, Mizukami M. A new finite element formulation for computational fluid dynamics: II. Beyond SUPG. Comput Methods Appl Mech Eng 1986;54:341–55.

[8] Fortin M. Old and new finite elements for incompressible flows. Int J Numer Methods Fluids 1981;1:347–64. [9] Verfürth R. Finite element approximation of incompressible Navier–Stokes equations with slip boundary condition II. Numer Math 1991;59:615–36. [10] Heywood JG, Rannacher R. 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 1982;19:275–311. [11] Heywood JG, Rannacher R. Finite element approximation of the nonstationary Navier–Stokes Problem. III. Smoothing property and higher order error estimates for spatial discretization. SIAM J Numer Anal 1988;25:489–512. [12] Bassi F, Crivellini A, Pietro DD, Rebay S. An implicit high-order discontinuous Galerkin method for steady and unsteady incompressible flows. Comput Fluids 2007;36:1529–46. [13] Shahbazi K, Fischer PF, Ethier CR. A high-order discontinuous Galerkin method for the unsteady incompressible Navier–Stokes equations. J Comput Phys 2007;222:391–407. [14] Ferrer E, Willden R. A high order discontinuous Galerkin finite element solver for the incompressible Navier–Stokes equations. Comput Fluids 2011;46:224–30. [15] Nguyen N, Peraire J, Cockburn B. An implicit high-order hybridizable discontinuous Galerkin method for the incompressible Navier–Stokes equations. J Comput Phys 2011;230:1147–70. [16] Rhebergen S, Cockburn B. A spacetime hybridizable discontinuous Galerkin method for incompressible flows on deforming domains. J Comput Phys 2012;231:4185–204. [17] Rhebergen S, Cockburn B, van der Vegt JJ. A spacetime discontinuous Galerkin method for the incompressible Navier–Stokes equations. J Comput Phys 2013;233:339–58. [18] Crivellini A, D’Alessandro V, Bassi F. High-order discontinuous Galerkin solutions of three-dimensional incompressible RANS equations. Comput Fluids 2013;81:122–33. [19] Klein B, Kummer F, Oberlack M. A SIMPLE based discontinuous Galerkin solver for steady incompressible flows. J Comput Phys 2013;237:235–50. [20] Bassi F, Crivellini A, Pietro DD, Rebay S. On a robust discontinuous Galerkin technique for the solution of compressible flow. J Comput Phys 2006;218:208–21. [21] Chorin A. A numerical method for solving incompressible viscous flow problems. J Comput Phys 1967;2:12–26.

M. Tavelli, M. Dumbser / Computers & Fluids 119 (2015) 235–249 [22] Chorin A. Numerical solution of the Navier–Stokes equations. Math Comput 1968;23:341–54. [23] Baumann C, Oden J. A discontinuous hp finite element method for convection– diffusion problems. Comput Methods Appl Mech Eng 1999;175(3–4):311–41. [24] Baumann C, Oden J. A discontinuous hp finite element method for the Euler and Navier–Stokes equations. Int J Numer Methods Fluids 1999;31(1):79–95. [25] Arnold D, Brezzi F, Cockburn B, Marini L. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J Numer Anal 2002;39:1749–79. [26] Tumolo G, Bonaventura L, Restelli M. A semi-implicit, semi-Lagrangian, padaptive discontinuous Galerkin method for the shallow water equations. J Comput Phys 2013;232:46–67. [27] Giraldo FX, Restelli M. High-order semi-implicit time-integrators for a triangular discontinuous Galerkin oceanic shallow water model. Int J Numer Methods Fluids 2010;63:1077–102. [28] Dolejsi V. Semi-implicit interior penalty discontinuous Galerkin methods for viscous compressible flows. Commun Comput Phys 2008;4:231–74. [29] Dolejsi V, Feistauer M. A semi-implicit discontinuous Galerkin finite element method for the numerical solution of inviscid compressible flow. J Comput Phys 2004;198:727–46. [30] Dolejsi V, Feistauer M, Hozman J. Analysis of semi-implicit DGFEM for nonlinear convection–diffusion problems on nonconforming meshes. Comput Methods Appl Mech Eng 2007;196:2813–27. [31] Liu YJ, Shu CW, Tadmor E, Zhang M. Central discontinuous Galerkin methods on overlapping cells with a non-oscillatory hierarchical reconstruction. SIAM J Numer Anal 2007;45:2442–67. [32] Liu YJ, Shu CW, Tadmor E, Zhang M. L2-stability analysis of the central discontinuous Galerkin method and a comparison between the central and regular discontinuous Galerkin methods. Math Model Numer Anal 2008;42:593–607. [33] Chung E, Engquist B. Optimal discontinuous Galerkin methods for wave propagation. SIAM J Numer Anal 2006;44:2131–58. [34] Chung ET, Lee CS. A staggered discontinuous Galerkin method for the convection–diffusion equation. J Numer Math 2012;20:1–31. [35] Dumbser M, Casulli V. A staggered semi-implicit spectral discontinuous Galerkin scheme for the shallow water equations. Appl Math Comput 2013;219(15):8057–77. [36] Tavelli M, Dumbser M. A high order semi-implicit discontinuous Galerkin method for the two dimensional shallow water equations on staggered unstructured meshes. Appl Math Comput 2014;234:623–44. [37] Tavelli M, Dumbser M. A staggered arbitrary high order semi-implicit discontinuous Galerkin method for the two dimensional incompressible Navier–Stokes equations. Appl Math Comput 2014;248:70–92. [38] Hirt CW, Nichols BD. Volume of fluid (VOF) method for dynamics of free boundaries. J Comput Phys 1981;39:201–25. [39] Casulli V, Cheng RT. Semi-implicit finite difference methods for threedimensional shallow water flow. Int J Numer Methods Fluids 1992;15:629–48. [40] Casulli V. A semi-implicit finite difference method for non-hydrostatic freesurface flows. Int J Numer Methods Fluids 1999;30:425–40. [41] Casulli V, Walters RA. An unstructured grid, three-dimensional model based on the shallow water equations. Int J Numer Methods Fluids 2000;32:331–48. [42] Casulli V. A high-resolution wetting and drying algorithm for free-surface hydrodynamics. Int J Numer Methods Fluids 2009;60:391–408. [43] Casulli V, Stelling GS. Semi-implicit subgrid modelling of three-dimensional free-surface flows. Int J Numer Methods Fluids 2011;67:441–9. [44] Casulli V. A semi-implicit numerical method for the free-surface Navier– Stokes equations. Int J Numer Methods Fluids 2014;74:605–22. [45] Klein R, Botta N, Schneider T, Munz C, Roller S, Meister A, et al. Asymptotic adaptive methods for multi-scale problems in fluid mechanics. J Eng Math 2001;39:261–343. [46] Roller S, Munz C. A low mach number scheme based on multi-scale asymptotics. Comput Vis Sci 2000;3:85–91. [47] Park J, Munz C-D. Multiple pressure variables methods for fluid flow at all mach numbers. Int J Numer Methods Fluids 2005;49(8):905–31. [48] van der Vegt JJW, van der Ven H. Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flows I. General formulation. J Comput Phys 2002;182:546–85. [49] van der Ven H, van der Vegt JJW. Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flows II. Efficient flux quadrature. Comput Methods Appl Mech Eng 2002;191:4747–80. [50] Klaij C, der Vegt JJWV, der Ven HV. Space-time discontinuous Galerkin method for the compressible Navier–Stokes equations. J Comput Phys 2006;217:589–611. [51] Demkowicz L, Oden J, Strouboulis T. Adaptive finite elements for flow problems with moving boundaries. Part I: Variational principles and a posteriori estimates. Comput Methods Appl Mech Eng 1984;46:217–51.

249

[52] Cockburn B, Shu CW. The Runge–Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems. J Comput Phys 1998;141:199–224. [53] Cockburn B, Shu CW. The local discontinuous Galerkin method for timedependent convection diffusion systems. SIAM J Numer Anal 1998;35:2440–63. [54] Cockburn B, Shu CW. Runge–Kutta discontinuous Galerkin methods for convection-dominated problems. J Sci Comput 2001;16:173–261. [55] Rusanov VV. Calculation of interaction of non-steady shock waves with obstacles. J Comput Math Phys USSR 1961;1:267–79. [56] Gassner G, Lörcher F, Munz CD. A contribution to the construction of diffusion fluxes for finite volume and discontinuous Galerkin schemes. J Comput Phys 2007;224:1049–63. [57] Bermudez A, Dervieux A, Desideri J, Vazquez M. Upwind schemes for the twodimensional shallow water equations with variable depth using unstructured meshes. Comput Methods Appl Mech Eng 1998;155:49–72. [58] Bermúdez A, Ferrín J, Saavedra L, Vázquez-Cendón M. A projection hybrid finite volume/element method for low-Mach number flows. J Comput Phys 2014;271:360–78. [59] Toro EF, Hidalgo A, Dumbser M. FORCE schemes on unstructured meshes I: conservative hyperbolic systems. J Comput Phys 2009;228:3368–89. [60] Dumbser M. Arbitrary high order PNPM schemes on unstructured meshes for the compressible Navier–Stokes equations. Comput Fluids 2010;39:60–76. [61] Castro MJ, Gallardo JM, Parés C. High-order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. applications to shallow-water systems. Math Comput 2006;75:1103–34. [62] Parés C. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM J Numer Anal 2006;44:300–21. [63] Maso GD, LeFloch PG, Murat F. Definition and weak stability of nonconservative products. J Math Pures Appl 1995;74:483–548. [64] Castro MJ, LeFloch PG, Muñoz-Ruiz ML, Parés C. Why many theories of shock waves are necessary: convergence error in formally path-consistent schemes. J Comput Phys 2008;227:8107–29. [65] Rhebergen S, Bokhove O, van der Vegt JJW. Discontinuous Galerkin finite element methods for hyperbolic nonconservative partial differential equations. J Comput Phys 2008;227:1887–922. [66] Dumbser M, Castro M, Parés C, Toro E. ADER schemes on unstructured meshes for non-conservative hyperbolic systems: applications to geophysical flows. Comput Fluids 2009;38:1731–48. [67] Dumbser M, Hidalgo A, Castro M, Parés C, Toro EF. FORCE schemes on unstructured meshes II: non-conservative hyperbolic systems. Comput Methods Appl Mech Eng 2010;199:625–47. [68] Bassi F, Rebay S. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations. J Comput Phys 1997;131:267–79. [69] Dumbser M, Balsara DS, Toro EF, Munz CD. A unified framework for the construction of one-step finite-volume and discontinuous Galerkin schemes. J Comput Phys 2008;227:8209–53. [70] Casulli V, Zanolli P. High resolution methods for multidimensional advection– diffusion problems in free-surface hydrodynamics. Ocean Model 2005;10:137–51. [71] Tavelli M, Dumbser M, Casulli V. High resolution methods for scalar transport problems in compliant systems of arteries. Appl Numer Math 2013;74:62–82. [72] Saad Y, Schultz M. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci Stat Comput 1986;7:856–69. [73] Dumbser M, Facchini M. A local space–time discontinuous Galerkin method for Boussinesq-type equations. Appl Math Comput, in press. http://dx.doi.org/ 10.1016/j.amc.2015.06.052. [74] Jiang G, Shu CW. On a cell entropy inequality for discontinuous Galerkin methods. Math Comput 1994;62:531–8. [75] Barth T, Charrier P. Energy stable flux formulas for the discontinuous Galerkin discretization of first-order nonlinear conservation laws. Tech Rep NAS-01001, NASA; 2001. [76] Womersley J. Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. J Physiol 1955;127:553–63. [77] Bell JB, Coletta P, Glaz HM. A second-order projection method for the incompressible Navier–Stokes equations. J Comput Phys 1989;85:257–83. [78] Ghia U, Ghia KN, Shin CT. High-re solutions for incompressible flow using Navier–Stokes equations and multigrid method. J Comput Phys 1982;48:387–411. [79] Qu L, Norberg C, Davidson L, Peng S, Wang F. Quantitative numerical analysis of flow past a circular cylinder at reynolds number between 50 and 200. J Fluids Struct 2013;39:347–70.