Semi-implicit discontinuous Galerkin method for the solution of the compressible Navier-Stokes equations

(1)

P. Wesseling, E. O˜nate and J. P´eriaux (Eds) c

°TU Delft, The Netherlands, 2006

SEMI-IMPLICIT DISCONTINUOUS GALERKIN METHOD

FOR THE SOLUTION OF THE COMPRESSIBLE

NAVIER-STOKES EQUATIONS

V´ıt Dolejˇs´ı and Jiˇr´ı Hozman

Charles University Prague, Faculty of Mathematics and Physics Sokolovsk´a 83, 186 75 Prague, Czech Republic

e-mail: dolejsi@karlin.mff.cuni.cz, jhozmi@volny.cz

Key words: compressible Navier-Stokes equations, discontinuous Galerkin method, backward difference formulae, semi-implicit scheme

Abstract. We deal with a numerical solution of the compressible Navier-Stokes equa-tions with the aid of higher order schemes. We employ a combination of the discontinuous Galerkin finite element method for the space semi-discretization and a backward difference formula for the time discretization. Moreover, using a linearization of inviscid as well as viscous fluxes and applying a suitable explicit extrapolation for nonlinear terms we have to solve only a linear algebraic problem at each time step. Then we obtain an efficient numerical scheme which is almost unconditionaly stable and has a higher degree of ap-proximation with respect to the space and time coordinates. Several numerical examples are presented.

1

1 INTRODUCTION

Our aim is to develop an efficient, robust and accurate numerical scheme for a sim-ulation of unsteady compressible flow. Very promising it seems to carry out the space discretization by the discontinuous Galerkin method (DGM), which is based on a piece-wise polynomial but discontinuous approximation. For a survey about DGM, see [7] or [9]. There are several variants of the DGM for the solution of problems containing diffusion terms. We prefer the discontinuous Galerkin finite elements (DGFE) method with a non-symmetric variant of stabilization and an interior and boundary penalty. This scheme is usually denoted as NIPG (nonsymetric interior penalty Galerkin) method. The advantage of this approach is that the corresponding diffusive form has a coercivity property for any positive value of a penalty coefficient σ, see the analysis in [1], [2], [12], [18], [19]. This is important for the case of the compressible Navier-Stokes equations where a numerical analysis is impossible.

1

(2)

It is possible to use a discontinuous approximation also for the time discretization [21], [22], but the most usual approach is an application of the method of lines. In this case, the Runge-Kutta methods are very popular for their simplicity and a high order of accuracy, see [3], [4], [5], [6], [8], [17]. Their drawback is a strong restriction to the choice of the time step. To avoid this disadvantage it is suitable to use an implicit time discretization. A full implicit scheme leads to a necessity to solve a nonlinear system of algebraic equations at each time step which is rather expensive. Therefore, we proposed in [14] a semi-implicit method for the simulation of inviscid compressible flow. This technique is based on a suitable linearization of the Euler fluxes. The linear terms are treated implicitly whereas the nonlinear ones explicitly which leads to a linear algebraic problem at each time step. In this paper we extend the approach of semi-implicit scheme to the viscous case where the diffusive terms are discretized by NIPG approach presented in [10]. So that this article is a natural combination of the explicit scheme for viscous flow from [10] with a semi-implicit scheme for inviscid flow from [14]. Moreover, we apply a backward difference formula (BDF) to the time discretization which gives a higher order approximation with respect to the time.

The contents of the rest of the paper is the following. In Section 2 we introduce the com-pressible Navier-Stokes equations with appropriate initial and boundary conditions and mention some properties of inviscid and viscous fluxes. Section 3 describes a discretization of the Navier-Stokes equations with the aid of DGFE method and a linearization of the inviscid and the viscous fluxes, which yields to a higher order semi-implicit scheme. Two numerical examples are presented in Section 4 and the concluding remarks are given in Section 5.

2 PROBLEM FORMULATION

2.1 Compressible flow problem

Let Ω ⊂ IR2 _{be a bounded domain and T > 0. We set Q}

T = Ω × (0, T ) and by ∂Ω

denote the boundary of Ω which consists of several disjoint parts. We distinguish inlet ΓI,

outlet ΓO and impermeable walls ΓW on ∂Ω. The system of the Navier-Stokes equations

describing a 2D motion viscous compressible flow can be written in the dimensionless form ∂w ∂t + 2 X s=1 ∂fs(w) ∂xs = 2 X s=1 ∂Rs(w, ∇w) ∂xs in QT, (1) where w = (w1, . . . , w4)T = (ρ, ρv1, ρv2, e)T (2)

is the so-called state vector,

fs(w) = (fs(1)(w), . . . , f (4) s (w))

T ₍₃₎

(3)

are the so-called inviscid (Euler) fluxes and Rs(w, ∇w) = (R(1)s (w, ∇w), . . . , R(4)s (w, ∇w))T (4) = Ã 0, τs1, τs2, 2 X k=1 τskvk+ γ Re P r ∂θ ∂xs !T , s = 1, 2

are the so-called viscous fluxes. We consider the Newtonian type of fluid, i. e., the viscous part of the stress tensor has the form

τsk = 1 Re "Ã ∂vs ∂xk + ∂vk ∂xs ! − 2 3div(v) δsk # , s, k = 1, 2. (5)

We use the following notation: ρ – density, p – pressure, e – total energy, v = (v1, v2) –

velocity, θ – temperature, γ – Poisson adiabatic constant, Re – Reynolds number, P r – Prandtl number.

In order to close the system we consider the following thermodynamical relations: the state equation for perfect gas and the relation for total energy

p = (γ − 1) (e − ρ|v|2/2), e = cVρθ + ρ|v|2/2, (6)

where cV is the specific heat at constant volume which is equal to one in the dimensionless

case. The system (1) – (6) is hyperbolic-parabolic type. It is equipped with the initial condition

w(x, 0) = w0(x), x ∈ Ω, (7)

and the following set of boundary conditions on appropriate parts of boundary: a) ρ = ρD, v = vD, 2 X k=1 Ã ₂ X l=1 τlknl ! vk+ γ Re P r ∂θ ∂n = 0 on ΓI, (8) b) 2 X k=1 τsknk= 0, s = 1, 2, ∂θ ∂n = 0 on ΓO, c) v = 0, ∂θ ∂n = 0 on ΓW,

where ρD and vD are given function and n = (n1, n2) is a unit outer normal to ∂Ω.

Another possibility is to replace the adiabatic boundary condition (8), c) by

c0) v = 0, θ = θD on ΓW. (9)

The problem to solve the compressible Navier-Stokes equations (1) – (4) with consti-tutive relations (5) – (6), equipped with the initial and boundary conditions (7) – (9) will be denoted by (CFP) (compressible flow problem).

In next two sections we present some properties of the inviscid and viscous fluxes fs, s = 1, 2 and Rs, s = 1, 2 given by (3) and (4), respectively. These properties are

(4)

2.2 Properties of the inviscid fluxes

Using relations (2) – (3), we express the Euler fluxes fs, s = 1, 2, in terms of the

variables w1, . . . , w4 in the form

fs(w) =            ws+1 ws+1w2 w1 + δs1(γ − 1) ³ w4− w2 2+w32 2w1 ´ ws+1w3 w1 + δs2(γ − 1) ³ w4− w2 2+w32 2w1 ´ ws+1 w1 ³ γw4− (γ − 1)w 2 2+w32 2w1 ´            , s = 1, 2. (10)

Obviously, fs, s = 1, 2, are homogeneous mappings of order one, i.e.,

fs(ξw) = ξfs(w), ξ ∈ IR, ξ 6= 0, s = 1, 2. (11)

Then it is easy to show (see [15], page 432) that

f_s(w) = As(w)w, s = 1, 2, (12)

where

As(w) ≡

Dfs(w)

Dw , s = 1, 2, (13)

are the Jacobi matrices of the mappings fs. Moreover we define a matrix

P(w, n) ≡

2

X

s=1

As(w)ns, (14)

where n = (n1, n2) ∈ IR2, n21+ n22 = 1, which plays a role in the definition of a numerical

flux and the choice of boundary conditions. 2.3 Properties of the viscous fluxes

The viscous terms Rs(w, ∇w) can be expressed in the form

Rs(w, ∇w) = 2 X k=1 Ksk(w) ∂w ∂xk , s = 1, 2, (15)

where Ksk are 4 × 4 matrices dependent on w, see, e.g., [16]. However, the numerical

experiments carried out in [10] indicate, that a linearization based on (15) gives unstable results. Therefore we introduced in [10] more sophisticated relations replacing (15). Let w= (w1, . . . , w4) ∈ IR4 and ϕ = (ϕ1, . . . , ϕ4) ∈ IR4, then we define the forms Ds(·, ·, ·, ·) :

(5)

(6)

D2,2(w) ≡       0 0 0 0 0 _{Re w}1 ₁ 0 0 0 0 4₃_{Re w}1 ₁ 0 0 w2 w1 1 Re w1 ³ 1 − γ P r ´ _w 3 w1 1 Re w1 ³ 4 3 − γ P r ´ _γ P r 1 Re w1       . (22)

It is possible to show that the forms Ds are consistent with viscous fluxes Rs in the

way

Ds(w, ∇w, w, ∇w) = Rs(w, ∇w) ∀w, s = 1, 2. (23)

Moreover, the forms Ds(w, ∇w, ϕ, ∇ϕ), s = 1, 2 are linear with respect ϕ and ∇ϕ and

they are independent of ∇ϕ1, whereas the form

Rs(w, ∇ϕ) = 2 X k=1 Ksk(w) ∂ϕ ∂xk , s = 1, 2, (24)

following from (15) depends on ∇ϕ1. This was the originality of the approach presented

in [10].

3 DISCRETIZATION

3.1 Broken Sobolev space

Let Th (h > 0) denote a triangulation of the closure Ω of the domain Ω into a finite

number of closed elements (triangles or quadrilaterals) K with mutually disjoint interiors. We set h = maxK∈Thdiam(K). Let I be a suitable index set such that Th = {Ki}i∈I.

If two elements Ki, Kj ∈ Th contain a nonempty open part of their faces, we put Γij =

Γji = ∂Ki∩ ∂Kj. For i ∈ I we set s(i) = {j ∈ I; Γij exists}. The boundary ∂Ω is formed

by a finite number of faces of elements Ki adjacent to ∂Ω. We denote all these boundary

faces by Sj, where j ∈ Ib is a suitable index set and put γ(i) = {j ∈ Ib; Sj is a face of Ki},

Γij = Sj for Ki ∈ Th such that Sj ⊂ ∂Ki, j ∈ Ib. For Ki not containing any boundary

face Sj we put γ(i) = ∅. Further, we define three disjoint subsets γI(i), γO(i) and γW(i)

corresponding to the boundary parts ΓI, ΓO and ΓW, respectively. Obviously, γ(i) =

γI(i) ∪ γO(i) ∪ γW(i), i ∈ I. Moreover, for simplicity we put γIO(i) ≡ γI(i) ∪ γO(i), i ∈ I

denoting the inflow/outflow parts of boundary and γD(i) ≡ γI(i) ∪ γW(i), i ∈ I denoting

the parts of boundary, where at least one Dirichlet boundary condition (8) is prescribed (inflow and solid wall). Finally, we put S(i) = s(i) ∪ γ(i) and nij = ((nij)₁, (nij)₂) is the

unit outer normal to ∂Ki on the face Γij.

Over the triangulation Th we define the broken Sobolev space

Hk(Ω, Th) = {v; v|K ∈ Hk(K) ∀ K ∈ Th}, (25)

where Hk_{(K) = W}k,2_{(K) denotes the (classical) Sobolev space on element K. For v ∈}

H1_{(Ω, T}

h), j ∈ s(i), i ∈ I we set

v|Γij = trace of v|Ki on Γij, (26)

(7)

denoting the traces of v on Γij = Γji, which are different in general. Moreover, [v]_Γ_ij = v|_Γ_ij− v|_Γ_ji (27) hvi_Γ ij = 1 2 ³ v|_Γ ij + v|Γji ´ (28) denotes the jump and mean value of function v over the edge Γij, respectively. Finally,

we put

hvi_Γ_ij = v|_Γ_ij for j ∈ γ(i), i ∈ I. (29) 3.2 Weak DGFE formulation

In the same way as in [10] we introduce the NIPG formulation of (CFP). For w, ϕ ∈ H2_{(Ω, T}

h)4 we define the forms

(w, ϕ) = X Ki∈Th Z Ki w· ϕ dx (30) (L2_{-scalar product),} ˜ ah(w, ϕ) = X Ki∈Th (Z Ki 2 X s=1 Rs(w, ∇w) · ∂ϕ ∂xs dx (31) − X j∈s(i) j<i Z Γij 2 X s=1 ³ hDs(w, ∇w, w, ∇w)i (nij)s· [ϕ] −hDs(w, ∇w, ϕ, ∇ϕ)i (nij)s· [w] ´ dS − X j∈γD(i) Z Γij 2 X s=1 ³ Ds(w, ∇w, w, ∇w) (nij)s· ϕ dS − Ds(w, ∇w, ϕ, ∇ϕ) (nij)s· (w − wB) ´ dSo (diffusion form), ˜ bh(w, ϕ) = X i∈I    X j∈S(i) Z Γij 2 X s=1 fs(w)ns· ϕ dS − Z Ki 2 X s=1 fs(w) · ∂ϕ ∂xs dx    (32)

(convective form), and

(8)

(interior and boundary penalty terms), where the penalty parameter σ is chosen as σ|Γij =

(|Γij| Re )−1. The state vector wB prescribed on ΓI ∪ ΓW is given by the boundary

conditions, in particular, for the case (8) a)–c) we have

wB = (ρ|ΓW, 0, 0, ρ|ΓWθ|ΓW) on ΓW, (34) wB = (ρD, ρD(vD)1, ρD(vD)2, ρ|ΓIθ|ΓI + 1 2ρD|vD| 2_{) on Γ} I,

and for the case (8) a)–b), (9) c’)

wB = (ρ|ΓW, 0, 0, ρ|ΓWθD) on ΓW, (35) wB = (ρD, ρD(vD)1, ρD(vD)2, ρ|ΓIθ|ΓI + 1 2ρD|vD| 2_{) on Γ} I,

where ρD, vD and θD are given functions from the boundary conditions (8)–(9) and ρ|Γ

and θ|Γ are the values of density and temperature extrapolated from interior of Ω on the

appropriate boundary part Γ, respectively.

Let w(t) denotes the function on Ω such that w(t) (x) = w(x, t), x ∈ Ω. Then with the aid of (30) – (33) the NIPG formulation for the Navier-Stokes equations reads

d

dt(w(t), ϕ) + ˜ah(w(t), ϕ) + ˜bh(w(t), ϕ) + Jh(w(t), ϕ) = 0, (36) w(t), ϕ ∈ H2(Ω, Th)4, t ∈ (0, T ).

3.3 Space semi-discretization

Now we shall introduce the space semi-discretization of (36). To evaluate the boundary integrals in (32) we use the approximation

Z Γij 2 X s=1 f_s(w(t)) (nij)_s· ϕ dS ≈ Z Γij H(w(t)|Γij, w(t)|Γji, nij) · ϕ dS, (37)

where H is a numerical flux, w(t)|Γij and w(t)|Γji are the values of w on Γij considered

from the interior and the exterior of Ki, respectively, and at time t. For details, see, e.g.

[15] or [16]. Then with the aid of (32) and (37) we define the form ¯ bh(w, ϕ) ≡ X i∈I    X j∈S(i) Z Γij H(w|Γij, w|Γji, nij) · ϕ dS − Z Ki 2 X s=1 fs(‘w) · ∂ϕ ∂xs dx   . (38)

(9)

where λ1, . . . , λ4 are the eigenvalues of P . We define the “positive” and “negative” part

of P by

P±(w, n) = T Λ±T−1, Λ± = diag (λ±₁, . . . , λ±₄). (40) Then the Vijayasundaram numerical flux reads

HV S(w1, w2, n) ≡ P+ µ_w 1+ w2 2 , n ¶ w1+ P− µ_w 1+ w2 2 , n ¶ w2. (41)

Similarly as the Vijayasundaram numerical flux (41) we can apply, e.g., the Roe’s type schemes [20] having also a form suitable for a linearization.

It is necessary to specify the meaning of w(t)|Γji for j ∈ γ(i). We use an approach

known from an inviscid flow simulation. On ΓI ∪ ΓO we prescribe mn components of w

on Γij and extrapolate mp components from the interior of Ki to Γij, where mn is the

number of negative eigenvalues of matrix P (w, n) given by (14). Thus, we define

w|Γji = LRP (w|Γij, wD, nij), (42)

where LRP (·, ·, ·) represents a solution of local Riemann problem considered on edge Γij

and wD is a given state vector (e.g. from far-field boundary conditions). For details, see

[13].

On ΓW, an impermeability condition

v· n = 0 (43)

is prescribed. Then in virtue of (37) we put

Z Γij H(w(t)|Γij, w(t)|Γji, nij) · ϕ dS := Z Γij FW(w(t), nij) · ϕ dS, j ∈ γW(i), (44) where FW(w, n) ≡ (0, pn1, pn2, 0)T. (45)

The pressure p is expressed in the form

p = (γ − 1) (w4− (w22+ w23)/(2w1), (46)

following from (6) and (2) and extrapolated on Γij from Ki and n = (n1, n2) = nij.

The approximate solution of (CFP) is sought in the space of discontinuous piecewise polynomial functions Sh defined by

Sh ≡ Sh× Sh× Sh× Sh, (47)

Sh ≡ {v; v|K ∈ Pp(K) ∀ K ∈ Th},

where p is a positive integer and Pp_{(K) denotes the space of all polynomials on K of}

degree at most p. Obviously, Sh ⊂ H1(Ω, T ) and therefore the forms (30), (31), (33), (38)

(10)

Definition 1 Function wh is a semidiscrete solution of (CFP), if a) wh ∈ C1([0, T ]; Sh), (48) b) Ã ∂wh(t) ∂t , ϕh ! + ˜ah(wh(t), ϕh) + ¯bh(wh(t), ϕh) + Jh(wh(t), ϕh) = 0 ∀ ϕ_h ∈ Sh, ∀ t ∈ (0, T ), c) wh(0) = w0h, where w0

h ∈ Sh denotes an Sh-approximation of the initial condition w0 from (7).

Here C1([0, T ]; Sh) is the space of continuously differentiable mappings of the interval

[0, T ] into Sh. The problem (48), a) – c) exhibits a system of ordinary differential equations

for wh(t) which has to be discretized by a suitable ODE method. In order to avoid the

drawbacks mentioned in Section 1 (time step restriction and nonlinearity of the discretized problem) we follow the approach presented in [14] for the inviscid flow simulation, where a semi-implicit time discretization was developed. Therefore we carry out a linearization of the nonlinear forms ˜ah and ¯bh. Firstly, in Section 3.4 we recall the linearization of ¯bh

presented in [14] and then in Section 3.5 we introduce the linearization of ˜ah.

3.4 Linearization of inviscid terms By (38), for wh, ϕh ∈ Sh we have ¯bh(wh, ϕh) = − X K∈Th Z K 2 X s=1 fs(wh(x)) · ∂ϕh(x) ∂xs dx | {z } =:˜σ1(wh,ϕh) (49) + X Ki∈Th X j∈S(i) Z Γij H(wh|Γij, wh|Γji, nij) · ϕhdS | {z } =:˜σ2(wh,ϕh) .

The individual forms ˜σ1(·, ·) and ˜σ2(·, ·) will be linearized separately. For ˜σ1, we employ

the property (12) of the Euler fluxes and for ¯wh, wh, ϕh ∈ Sh define a form

σ1( ¯wh, wh, ϕh) ≡ X K∈Th Z K 2 X s=1 As( ¯wh(x))wh(x) · ∂ϕh(x) ∂xs dx. (50) Obviously, ˜ σ1(wh, ϕh) = σ1(wh, wh, ϕh) ∀wh, ϕh ∈ Sh. (51)

The linearization of the term ˜σ2 can be carried out in a simple way, when H in (49) is

(11)

we put ˜ σ2( ¯wh, wh, ϕh) ≡ X Ki∈Th X j∈S(i) Z Γij h P+(h ¯whiij, nij) wh|Γij (52) + P−_{(h ¯}_w hiij, nij) wh|Γji i · ϕhdS, where h ¯wk hiij is given by (28) – (29).

It is necessary to pay a special attention to wh|Γji for Γij ⊂ ∂Ω. For j ∈ γW(i) we use

approximation (44). The vector FW given by (45) is a nonlinear function of w and its

linearization is given with the aid of the Taylor expansion as

FW(wh, n) ≈ ˜FW( ¯wh, wh, n) ≡ FW( ¯wh, n) + DFW( ¯wh, n) (wh− ¯wh) = DFW( ¯wh, n)wh, (53) where DFW(w, n) ≡ (γ − 1)      0 0 0 0 (v2 1 + v22)n1/2 −v1n1 −v2n1 n1 (v2 1 + v22)n2/2 −v1n2 −v2n2 n2 0 0 0 0      (54)

is obtained by the differentiation of function FW given by (45) with respect to w =

(w1, . . . , w4). Here n = (n1, n2), vj = wj+1/w1, j = 1, 2. The last equality in (53) follows

from (12) which implies FW( ¯wh, n) = DFW( ¯wh, n) ¯wh.

Then we put σ2( ¯wh, wh, ϕh) (55) ≡ X Ki∈Th X j∈s(i) Z Γij h P+(h ¯whiij, nij) wh|Γij+ P − (hwhiij, nij) wh|Γji i · ϕ_hdS + X Ki∈Th X j∈γIO(i) Z Γij h P+(h ¯whiij, nij) ¯wh|Γij + P − (h ¯whiij, nij) ¯wh|Γji i · ϕhdS + X Ki∈Th X j∈γW(i) Z Γij ˜ FW( ¯wh, wh, nij) · ϕ dS.

In (55) we use an approximation of type P±(·, ·) ¯w for j ∈ γIO(i) since this approach is

more suitable for implementation. Obviously, ˜

σ2(wh, ϕh) = σ2(wh, wh, ϕh) ∀wh, ϕh ∈ Sh. (56)

Finally, we define the form

bh( ¯wh, wh, ϕh) ≡ −σ1( ¯wh, wh, ϕh) + σ2( ¯wh, wh, ϕh), (57)

where σ1 and σ2 are given by (50) and (55), respectively. The form bh is linear with

respect to the second and third variable and using (51) and (56) consistent with ¯bh by

¯

(12)

3.5 Linearization of viscous terms

In virtue of (16) we can rewrite the diffusive form ˜ah(·, ·) given by (31) in an equivalent

form ˜ ah(w, ϕ) (59) = X Ki∈Th (Z Ki 2 X s=1 Ã Ds,0(w, ∇w) w + 2 X k=1 Ds,k(w) ∂w ∂xk ! · ∂ϕ ∂xs dx − X j∈s(i) j<i Z Γij 2 X s=1 Ã * Ds,0(w, ∇w) w + 2 X k=1 Ds,k(w) ∂w ∂xk + (nij)s· [ϕ] − * Ds,0(w, ∇w) ϕ + 2 X k=1 Ds,k(w) ∂ϕ ∂xk + (nij)s· [w] ! dS − X j∈γD(i) Z Γij 2 X s=1 µ³ Ds,0(w, ∇w) w + 2 X k=1 Ds,k(w) ∂w ∂xk ´ (nij)s· ϕ −³Ds,0(w, ∇w) ϕ + 2 X k=1 Ds,k(w) ∂ϕ ∂xk ´ (nij)s· (w − wB) ¶ dS ) , where Ds,k, s = 1, 2, k = 0, 1, 2 are given by (17) – (22). Then based on (59), for

¯

wh, wh, ϕh ∈ Sh we define the form

ah( ¯wh, wh, ϕh) (60) ≡ X Ki∈Th (Z Ki 2 X s=1 Ã Ds,0( ¯wh, ∇ ¯wh) wh+ 2 X k=1 Ds,k( ¯wh) ∂wh ∂xk ! · ∂ϕ ∂xs dx − X j∈s(i) j<i Z Γij 2 X s=1 ³* Ds,0( ¯wh, ∇ ¯wh) wh+ 2 X k=1 Ds,k( ¯wh) ∂wh ∂xk + (nij)s· [ϕ] − * Ds,0( ¯wh, ∇ ¯wh) ϕ + 2 X k=1 Ds,k( ¯wh) ∂ϕ ∂xk + (nij)s· [wh] ´ dS − X j∈γD(i) Z Γij 2 X s=1 µ³ Ds,0( ¯wh, ∇ ¯wh) wh+ 2 X k=1 Ds,k( ¯wh) ∂wh ∂xk ´ (nij)s· ϕ −³Ds,0( ¯wh, ∇ ¯wh) ϕ + 2 X k=1 Ds,k( ¯wh) ∂ϕ ∂xk ´ (nij)s· (wh− wB) ¶ dS ) , which is linear with respect to its second and third components. Moreover it is consistent with ˜ah(·, ·) by

˜

(13)

3.6 Full space-time discretization

The main idea of the semi-implicit discretization is to treat the linear parts of forms ah and bh (represented by their second arguments) implicitly and their nonlinear parts

(represented by their first arguments) explicitly. In order to obtain a sufficiently accurate approximation with respect to the time coordinate we use the so-called backward difference formula (BDF) for the solution ODE problem (48), a) – c). Moreover, for the nonlinear parts of ah(·, ·, ·) and bh(·, ·, ·) we employ a suitable explicit higher order extrapolation

which preserve a given order of accuracy and does not destroy the linearity of the problem at each time level.

Let 0 = t0 < t1 < . . . < tr = T be a partition of the interval (0, T ) and τk ≡

tk+1− tk, k = 0, 1, . . . , r − 1.

Definition 2 We define the approximate solution of (CFP) as functions wk

h, k = 1, . . . , r,

satisfying the conditions

a) wk+1_h ∈ Sh, (62) b) 1 τk Ã _n X l=0 αlwk+1−lh , ϕh ! + ah Ã _n X l=1 βlwk+1−lh , w k+1 h , ϕh ! +bh Ã _n X l=1 βlwk+1−lh , wk+1h , ϕh ! + Jh ³ wk+1_h , ϕh ´ = 0 ∀ ϕh ∈ Sh, k = n − 1, . . . , r − 1, c) w0h is Sh approximation of w0, d) wl_h ∈ Sh, l = 1, . . . , n − 1 are given

by a suitable one-step method,

where n ≥ 1 is the degree of the BDF scheme, the coefficients αl, l = 0, . . . , n and

βl, l = 1, . . . , n depend on time steps τk−l, l = 0, . . . , n.

The relations for the coefficients αl, l = 0, . . . , n and βl, l = 1, . . . , n were derived in

[11] for n = 1, 2, 3 and their values for constant time step τk = τ, k = 1, . . . , r are given

in Table 1.

n αl, l= 0, . . . , n βl, l= 1, . . . , n

1 1, −1 1

2 3₂, −2, 1₂ 2, −1

3 11₆ , −3, 3₂, −1₃ 3, −3, 1

Table 1: Values of coefficients αl and βl for constant time step

(14)

Based on numerical experiments carried out for a scalar nonlinear convection-diffusion equation in [11] we suppose that the orders of convergence of scheme (62), a) – d) are

O(hp¯+ τn), p =¯ ( p for p even p + 1 for p odd (63) in L∞_{((0, T ); L}2_{(Ω))-norm and} O(hp+ τn) (64)

in L∞_{((0, T ); H}1_{(Ω)) for a sufficiently regular exact (weak) solution.}

4 NUMERICAL EXAMPLES

In this section we present two numerical examples. The first one is a basic example of steady viscous flow and the second one is a classical benchmark for unsteady inviscid flow simulation.

Simulation of viscous flow We consider the laminar flow on the adiabatic flat plate characterised by a freestream Mach number M = 0.3 and by a Reynolds number Re = 104_.

The computation has been performed with the aid of piecewise linear approximation on a unstructured grid having 3781 elements which was adaptive refined along the flat. Figure 1 shows the comparison of the computed skin friction coefficient cf with the “theoretical”

one which is given by the well-known Blasius formula for the cf distribution along a

flat plate for incompressible flow. The right figure represents the distribution plotted in logarithmic scale. The computed results show good agreement with the Blasius solution. Simulation of unsteady inviscid flow We consider a flow through the forward facing step proposed in [24]. The geometry of the problem is shown at Figure 2. The vertical straight parts of boundary located at x = 0 and x = 3 exhibits inflow/outflow part of ∂Ω and other parts are impermeable walls. The initial condition is constant with the values

ρ = 1.4, v = (3, 0), p = 1 (65)

and the values of inflow/outflow boundary conditions are taken from (65). Therefore the inlet Mach number is M = 3. The simulation was performed until t = 3.5 with the aid of the scheme (62) with ah ≡ 0 and Jh ≡ 0 since the flow is inviscid. We carried out

three simulations on different grids with different degree of polynomial approximation. Table 2 shows the degree of polynomial approximation, number of triangles (#Th) of the

appropriate mesh and the total number of degree of freedom(#dof ).

Figure 2 (at top) shows the coarser mesh used for piecewise cubic approximation. The three step BDF formulae (m = 3 in (62)) was used for all three simulations. Figure 2 shows a comparison of obtained results with P1, P2 and P3 approximations, namely

isolines of Mach number at t = 3.5. We observe quantitatively similar results, but the P3 approximation given slightly sharper capturing of shock waves. Figures 3 and 4 show

(15)

0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.2 0.4 0.6 0.8 1 ’Blasius’ ’P_1’ 0.001 0.01 0.1 0.01 0.1 1 ’Blasius’ ’P_1’

Figure 1: Distribution of the computed skin friction coefficient (P 1) with the “theoretical” one (Blasius). The left and right figures are plotted in linear and logarithmic scaling, respectively

degree #Th #dof

1 9 932 119 184 2 2 914 69 936 3 1 033 41 320

Table 2: Polynomial degrees, the numbers of triangles and the numbers of degree of freedoms

5 CONCLUSION

We carried out a numerical solution of the compressible Navier-Stokes equations by a combination of DGFE method and BDF. This scheme is theoretically unconditionally stable, have a high order of approximation with respect to space and time coordinates and at each time step we solve a linear algebraic problem. Preliminary numerical examples give promising results.

REFERENCES

[1] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified analysis of discon-tinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal., 39(5):1749– 1779, 2002.

(16)

[3] F. Bassi and S. Rebay. A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations. J. Comput. Phys., 131:267–279, 1997.

[4] F. Bassi and S. Rebay. High-order accurate discontinuous finite element solution of the 2D Euler equations. J. Comput. Phys., 138:251–285, 1997.

[5] F. Bassi and S. Rebay. A high order discontinuous Galerkin method for compressible turbulent flow. In B. Cockburn, G. E. Karniadakis, and C.-W. Shu, editors, Discon-tinuous Galerkin Method: Theory, Computations and Applications, Lecture Notes in Computational Science and Engineering 11, pages 113–123. Springer-Verlag, 2000. [6] C. E. Baumann and J. T. Oden. A discontinuous hp finite element method for the

Euler and Navier-Stokes equations. Int. J. Numer. Methods Fluids, 31(1):79–95, 1999.

[7] B. Cockburn. Discontinuous Galerkin methods for convection dominated problems. In T. J. Barth and H. Deconinck, editors, High–Order Methods for Computational Physics, Lecture Notes in Computational Science and Engineering 9, pages 69–224. Springer, Berlin, 1999.

[8] B. Cockburn, S. Hou, and C. W. Shu. TVB Runge-Kutta local projection discontin-uous Galerkin finite element for conservation laws IV: The multi-dimensional case. Math. Comp., 54:545–581, 1990.

[9] B. Cockburn, G. E. Karniadakis, and C.-W. Shu (eds.). Discontinuous Galerkin methods. In Lecture Notes in Computational Science and Engineering 11. Springer, Berlin, 2000.

[10] V. Dolejˇs´ı. On the discontinuous Galerkin method for the numerical solution of the Navier–Stokes equations. Int. J. Numer. Methods Fluids, 45:1083–1106, 2004. [11] V. Dolejˇs´ı. Higher order semi-implicit discontinuous Galerkin finite element schemes

for compressible flow simulation. In Software and Algorithms of Numerical Mathe-matics, (submitted) 2005.

[12] V. Dolejˇs´ı, M. Feistauer, and V. Sobot´ıkov´a. A discontinuous Galerkin method for nonlinear convection–diffusion problems. Comput. Methods Appl. Mech. Engrg., 194:2709–2733, 2005.

(17)

[14] V. Dolejˇs´ı and M. Feistauer. Semi-implicit discontinuous Galerkin finite element method for the numerical solution of inviscid compressible flow. J. Comput. Phys., 198(2):727–746, 2004.

[15] M. Feistauer. Mathematical Methods in Fluid Dynamics. Longman Scientific & Technical, Harlow, 1993.

[16] M. Feistauer, J. Felcman, and I. Straˇskraba. Mathematical and Computational Meth-ods for Compressible Flow. Oxford University Press, Oxford, 2003.

[17] R. Hartmann and P. Houston. Adaptive discontinuous Galerkin finite element meth-ods for the compressible Euler equations. J. Comput. Phys., 183(2):508–532, 2002. [18] J. T. Oden, I. Babuˇska, and C. E. Baumann. A discontinuous hp finite element

method for diffusion problems. J. Comput. Phys, 146:491–519, 1998.

[19] B. Rivi`ere and M. F. Wheeler. A discontinuous Galerkin method applied to nonlin-ear parabolic equations. In B. Cockburn, editor, Discontinuous Galerkin methods. Theory, computation and applications., volume 11 of Lect. Notes Comput. Sci. Eng., pages 231–244. Berlin, Springer, 2000.

[20] P. L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes. J. Comput. Phys., 43(2):357–372, 1981.

[21] J. J. W. van der Vegt and H. van der Ven. Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flows. I: General formulation. J. Comput. Phys., 182(2):546–585, 2002.

[22] H. van der Ven and J. J. W. van der Vegt. Space-time discontinuous Galerkin finite element method with dynamic grid motion for inviscid compressible flows II. efficient flux quadrature. Comput. Methods Appl. Mech. Engrg., 191:4747–4780, 2002.

[23] G. Vijayasundaram. Transonic flow simulation using upstream centered scheme of Godunov type in finite elements. J. Comput. Phys., 63:416–433, 1986.

(18)

0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3

mesh for P3-approximation with #Th = 1 033

P1-approximation on #Th = 9 932

(19)

t= 0.2

t= 0.5

t= 1.0

t= 1.5

(20)

t= 2.0

t= 2.5

t= 3.0

t= 3.5