A new technique for the numerical solution of the compressible Euler - - PowerPoint PPT Presentation

A new technique for the numerical solution of the compressible Euler equations with arbitrary Mach numbers ∗ Miloslav Feistauer and V´ aclav Kuˇ cera Charles University Prague Faculty of Mathematics and Physics

∗Presented at conference HYP06, Lyon 17–21 July 2006.

Our goal: to develop a sufficiently accurate and robust method allowing the solution of problems with a wide range

• f Mach numbers

Inviscid flow 2D flow:

Basic equations: continuity eq., Euler equtions of motion, en- ergy equation = nonlinear hyperbolic system of conservation laws ∂w ∂t +

2

• s=1

∂fs(w) ∂xs = 0 in QT = Ω × (0, T), (1)

w = (w1, . . . , w4)T = (ρ, ρv1, ρv2, E)T

(2)

fi(w)

= (fi1, . . . , fim)T (3) = (ρvi, ρv1vi + δ1i p, ρv2vi + δ2i p, (E + p)vi)T p = (γ − 1) (E − ρ|v|2/2). (4)

Notation:

ρ – density, p – pressure, E – total energy,

v = (v1, v2) – velocity, δsk – Kronecker symbol, γ > 1 – Pois-

son adiabatic constant

Initial condition:

w(x, 0) = w0(x),

x ∈ Ω, (5)

Boundary conditions chosen in such a way that the

problem is linearly well–posed

Important property of fluxes:

As(w) = Dfs(w)

Dw , (6) – the Jacobi matrices of the mappings fs

fs, – homogeneous mappings of order one =

⇒

fs(w) = As(w)w

(7)

Discontinuous Galerkin (DG) space discretization

for N = 2

Mesh:

Ωh = polygonal approximation of Ω Th = mesh in Ωh consisting of triangles or quadrilaterals Ki ∈ Th, i ∈ I (= suitable index set) Γij = common edge between two neighbouring elements Ki and Kj or face ⊂ Ωh Introduce index sets s(i) so that j ∈ s(i) = ⇒ Γij ⊂ Ωh γ(i) so that i ∈ γ(i) = ⇒ Γij ⊂ ∂Ωh S(i) = s(i) ∪ γ(i). Then ∂Ki =

• j∈S(i)

Γij, ∂Ki ∩ ∂Ωh =

• j∈γ(i)

Γij. (8)

nij = ((nij)1, (nij)2) = unit outer normal to ∂Ki on the side

Γij

K Γ n

i ij ij j

K

Neighbouring elements Ki, Kj – mesh can be nonconforming

Space of approximate solutions

Discontinuous piecewise polynomial functions: Sh = [Sh]4, (9) Sh ≡ Sr,−1(Ω, Th) = {v; v|K ∈ Pr(K) ∀ K ∈ Th}, r ≥ 0 – integer and Pr(K) denotes the space of all polyno- mials on K of degree ≤ r.

Derivation of the discrete problem:

w – exact sufficiently regular solution

multiply (1) by any ϕ ∈ Sh, integrate over any Ki, i ∈ I, apply Green’s theorem, sum over all Ki, i ∈ I

= ⇒ ∂ ∂t

• Ki∈Th
• Ki

w · ϕ dx

(10) −

• Ki∈Th
• Ki

2

• s=1

fs(w) · ∂ϕ

∂xs dx − +

• Ki∈Th
• j∈S(i)
• Γij

2

• s=1

fs(w) · ϕ ns dS = 0.

Numerical flux

Approximation of fluxes through interfaces Γij – with the aid of a numerical flux H(α, β, n): let wh, ϕh ∈ Sh

• Γij

2

• s=1

fs(wh) ns · ϕh dS

(11) ≈

• Γij

H(wh|Γij, wh|Γji, nij) · ϕh dS wh|Γij = the value of wh on Γij considered from the interior

• f Ki, wh|Γji = the value of wh on Γij considered from the

exterior of Ki

Properties of the numerical flux:

– (locally) Lipschitz continuous – consistent: H(α, α, n) = N

s=1 fs(α)ns

– conservative: H(α, β, n) = −H(β, α, −n)

Definition of forms

For wh, ϕh ∈ Sh we define (wh, ϕh)h =

• Ωh

wh · ϕh dx,

(12) Bh(wh, ϕh) (13) = −

• K∈Th
• K

2

• s=1

fs(wh) · ∂ϕh

∂xs dx +

• Ki∈Th
• j∈S(i)
• Γij

H(wh|Γij, wh|Γji, nij) · ϕhdS.

Approximate DGFE semidiscrete solution: a)

wh ∈ C1([0, T], Sh),

b) d dt (wh(t), ϕh) + Bh(wh(t), ϕh) = 0 (14) ∀ ϕh ∈ Sh ∀ t ∈ (0, T), c)

wh(0) = Πhw0 = L2−projection of w0 on Sh

(14) ≡ large system of ordinary differential equations – usually solved by explicit Runge-Kutta methods – work well for fast flow - conditionally stable!!! Slow flows – with very low Mach numbers – represent big

• bstacle for all existing methods: FD, FV, FE

− → general opinion that this type of flow can be solved

• nly on the basis of modified governing equations

using physical variables and derived on the basis of a multiscale analysis

• R. Klein, C.-D. Munz, A. Meister,...

Question: Is it possible to develop a method using conser- vation form of the Euler equations and robust with respect to the Mach number?

Successful technique allowing the solution of all Mach numbers flows: DGFEM

& Ingredients: 1) Semi-implicit time discretization (M.F.& V. Dolejˇ s ´ ı, JCP 2004): For each k ≥ 0 find wk+1

h

such that a)

wk+1

h

∈ Sh, (15) b)

 wk+1

h

− wk

h

τ , ϕh

 

h

+ bh(wk

h, wk+1 h

, ϕh) = 0, ∀ ϕh ∈ Sh, k = 0, 1, . . . , c)

w0

h = Πhw0.

bh(wk

h, wk+1 h

, ϕh) - linear with respect to ϕh and wk+1

h

Definition of bh: Partial linearization of Bh

We have Bh(wk+1

h

, ϕh) (16) = −

• K∈Th
• K

2

• s=1

fs(wk+1

h

(x)) · ∂ϕh(x) ∂xs dx

• =:˜

σ1

+

• Ki∈Th
• j∈S(i)
• Γij

H(wk+1

h

|Γij, wk+1

h

|Γji, nij) · ϕhdS

• =:˜

σ2

.

Linearization of ˜ σ1:

Homogeneity of fs ⇒

2

• s=1

fs(w)ns =

2

• s=1

As(w)wns.

(17) − → ˜ σ1 ≈ σ1 (18) =

• K∈Th
• K

2

• s=1

As(wk

h(x))wk+1 h

(x) · ∂ϕh(x) ∂xs dx.

Linearization of ˜ σ2 =

• Ki∈Th
• j∈S(i)
• Γij H(wk+1

h

|Γij, wk+1

h

|Γji, nij) · ϕhdS: Assume that

H = Vijayasundaram numerical flux: The matrix P (w, n) =

2

s=1 As(w)ns is diagonalizable:

P (w, n) = T DT −1, D = diag (λ1, . . . , λ4),

(19) λ1, . . . , λ4 = the eigenvalues of P . “Positive” and “negative” part of P :

P ±(w, n) = T D±T −1, D± = diag (λ±

1 , . . . , λ± 4 ).

(20) Vijayasundaram numerical flux:

HV (w1, w2, n)

(21) = P +

w1 + w2

2 , n

• w1 + P −

w1 + w2

2 , n

• w2.

= ⇒ ˜ σ2 ≈

• Ki∈Th
• j∈S(i)
• Γij
• P +

wk

hij, nij

• wk+1

h

|Γij (22) + P − wk

hij, nij

• wk+1

h

|Γji

• · ϕhdS,

where wk

hij ≡

1 2

• wk

h|Γij + wk h|Γji

• .

(23)

= ⇒ Linearized convection form: bh(wk

h, wk+1 h

, ϕh) (24) = −

• K∈Th
• K

2

• s=1

As(wk

h(x))wk+1 h

(x) · ∂ϕh(x) ∂xs dx +

• Ki∈Th
• j∈S(i)
• Γij
• P +

wk

hij, nij

• wk+1

h

|Γij + P − wk

hij, nij

• wk+1

h

|Γji

• · ϕhdS,

which is linear with respect to wk+1

h

and ϕh.

Semi-implicit linearized numerical scheme: for each

k ≥ 0 find wk+1

h

such that a)

wk+1

h

∈ Sh, (25) b)

• wk+1

h

, ϕh

• h + τkbh(wk

h, wk+1 h

, ϕh) =

• wk

h, ϕh

• h

∀ ϕh ∈ Sh, k = 0, 1, . . . , c)

w0

h = Πhw0.

2) The use of isoparametric elements in order to obtain an accurate solution near curved boundaries 3) Special treatment of boundary conditions transparent for accoustic waves, based on the use of characteristic vari- ables 4) Shock capturing for avoiding the Gibbs phenomenon manifested by spurious overshoots and undershoots in com- puted quantities near discontinuities (shock waves, contact discontinuities) Artificial viscosity proposed by Jaffr´ e, Johnson and Szepessy

• introduces some amount of artificial viscosity everywhere
• nonphysical entropy production

We introduce stabilization terms combining ideas of Jaffr´ e, Johnson and Szepessy and M.F., Dolejˇ s ´ ı and Schwab.

a)Define the discontinuity indicator gk(i) proposed by M.F., Dolejˇ s ´ ı and Schwab: gk(i) =

• ∂Ki

[ρk

h]2 dS/(hKi|Ki|3/4),

Ki ∈ Th. (26) [u]|Γij = uΓij − uΓji = the jump on Γij of a function u ∈ Sh. b)Define the discrete indicator Gk(i) = 0 if gk(i) < 1, Gk(i) = 1 if gk(i) ≥ 1, Ki ∈ Th. (27) c)To the left-hand side of (15), b) we add the artificial viscosity form βh(wk

h, wk+1 h

, ϕ) = ν1

• i∈I

hKiGk(i)

• Ki

∇wk+1

h

· ∇ϕ dx (28) with ν1 ≈ 1. d)Augment the left-hand side of (15), b) by adding the form Jh(wk

h, wk+1 h

, ϕ) = ν2

• i∈I
• j∈s(i)

1 2(Gk(i)+Gk(j))

• Γij

[wk+1

h

]·[ϕ] dS, (29)

where ν2 ≈ 1. Resulting scheme: a)

wk+1

h

∈ Sh, (30) b)

 wk+1

h

− wk

h

τk , ϕh

 

h

+ bh(wk

h, wk+1 h

, ϕh) +βh(wk

h, wk+1 h

, ϕh) + Jh(wk

h, wk+1 h

, ϕh) = 0, ∀ ϕh ∈ Sh, k = 0, 1, . . . , c)

w0

h = Πhw0.

This method successfully overcomes problems with the Gibbs phenomenon in the context of the semi-implicit scheme. Important: Gk(i) vanishes in regions where the solution is regular. = ⇒ The scheme does not produce any nonphysical entropy in these regions.

Numerical examples

In order to show the robustness of the described technique with respect to the Mach number, we present computa- tional results obtained for 3 types of stationary compress- ible flow (obtained by the time stabilization for ”t → ∞”). Piecewise quadratic elements (r = 2) are used. 1) Low Mach number flow past a circular cylinder far field velocity parallel to the axis x1 far field Mach number M∞ = 10−4

a) Coarse mesh – 16 vertices on the cylinder, 361 elements in the computational domain Velocity isolines for the approximate solution of compressible flow (left) and for the exact solution of incompressible flow (right)

Velocity distribution along the cylinder .

–1 1 2 3 4

b) Fine mesh – 364 vertices on the cylinder, 8790 elements in the computational domain Velocity isolines for the approximate solution of compressible flow (left) and for the exact solution of incompressible flow (right)

Velocity distribution along the cylinder (full line – compressible flow, dotted line – icompressible flow)

0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 –1 1 2 3 4

Stability of the scheme:

the method is practically unconditionally stable: Choice of CFL: start with CFL = 38, gradually increased up to 2000 ρmax − ρmin = 2.3 · 10−8, maxK∈Th|∇ρh|K| < 1.99 · 10−6 = ⇒ the density is practically constant = ⇒ the flow behaves as incompressible

2) Low Mach number rotational flow past a half-cylinder (recommended by P. Roe) exact solution of incompressible flow by L. Fraenkel, J. Fluid Mech. 1961: v∞ = (x2, 0) Rotational incompressible flow - streamlines of the exact solution Rotational compressible flow with M∞ = 0.0001 - streamlines of the approximate solution

3) Transonic flow Flow through the GAMM channel with a 10% circular bump and the inlet Mach number equal to 0.67

0.2 0.4 0.6 0.8 1

Transonic flow through the GAMM channel, Mach number isolines

Density distribution on the lower wall Zierep singularity

Important criterion of the quality of the method: the dis- tribution of entropy – produced on discontinuities only

0.2 0.4 0.6 0.8 1

Entropy distribution

Conclusion

Obtained results show that DGFEM is a promising robust and accurate technique for the numerical solution of singularly perturbed problems with dominating convection, conservation laws and com- pressible flow. Suitable space-time discretization, the treatment of bound- ary conditions and stabilization lead to a method, which allows the solution of compressible flow with a wide range

• f Mach numbers, using the conservative variables without

any modification of the governing equations.

Further goals: – extension to viscous compressible flow with a wide range of Mach numbers – development of the theoretical justification of the proposed method – application of the DGFEM to the numerical simulation of com- plex technically relevant problems – extension to 3D problems – application of higher degree elements (r > 2) – developement of a posteriori error estimates and space-time adaptivity (for nonstationary flow) – extension to flows in time dependent domains and aeroelastic problems

Thank you very much for your attention