WENO Shock-Fitted Solution to 1-D Euler Equations with Reaction - - PowerPoint PPT Presentation

WENO Shock-Fitted Solution to 1-D Euler Equations with Reaction Andrew K. Henrick, Tariq D. Aslam, Joseph M. Powers University of Notre Dame Los Alamos National Laboratory Tenth SIAM International Conference on Numerical Combustion Sedona, Arizona 10 May 2004

Introduction

• Euler equations with reaction.
• Simple two species irreversible exothermic reaction.
• High order numerical method to resolve behavior.
• Shock fitting to maintain order in presence of

discontinuity.

• ZND structure for f=1.8 and 1.6.

Review

• Fickett and Davis, Detonation: Theory and

Experiment, 1979.

• LeVeque, Numerical Methods for Cons. Laws, 1992
• Osher and Fedkiw, L.S.M. Dyn. Imp. Surfaces, 2003
• Jiang and Shu, J. Comp. Phys., 1996
• Short and Sharpe - 1D linear stability results
• Bdzil - Shock Change Equation

Numerical Difficulties

• PDE versus conservation at the shock.
• Shock capturing - oscillations.
• Shock tracking - evolution of shock location on fixed

grid.

• Shock fitting - numerical grid aligned with shock locus.

Governing Equations in 1-D

∂ρ ∂t + ∂ ∂x (ρu) = 0 ∂ ∂t (ρu) + ∂ ∂x

• ρu2 + p
• = 0

∂ ∂t

• ρ
• e + u2

2

• + ∂

∂x

• u(ρ
• e + u2

2

• + p)
• = 0

∂ ∂t (ρλ) + ∂ ∂x (ρλu) = kρ(1 − λ) exp −Eρ p

• p = ρRT

e =

1 γ−1p/ρ − λq

Jump Conditions

νι V(1) V(2) γ(1) γ(2) σ(1) Σ σ(2)

Governing PDE’s do not hold over a discontinuity - must appeal to integral formulation.

d dt Z

M− V

Fd − V = Z

M− V

Bd − V | {z }

volume forces

+ Z

MS

TinidS | {z }

surface forces

F(Di − vi) + Ti = 0 ρ(D − u) = 0 ρu(D − u) − p = 0 ρ „ e + u2 2 « (D − u) − up = 0 ρλ(D − u) = 0

This gives ρ(D), u(D), p(D), λ(D) given the quiescent state.

Wave Frame Transformation

¯ x(x, t) = x − s(t) −∞ < x < +∞ ¯ t(x, t) = t

all x and t,

9 = ;

∂f ∂t

= ∂ ¯

f ∂¯ x ∂¯ x ∂t + ∂ ¯ f ∂¯ t ∂¯ t ∂t

= ∂ ¯

f ∂¯ x(− ˙

s(t)) + ∂ ¯

f ∂¯ t

= −D ∂ ¯

f ∂¯ x + ∂ ¯ f ∂¯ t ∂f ∂x

= ∂ ¯

f ∂¯ x ∂¯ x ∂x + ∂ ¯ f ∂¯ t ∂¯ t ∂x

= ∂ ¯

f ∂¯ x

∂ρ ∂¯ t + ∂ ∂¯ x (ρu − Dρ) = 0 ∂ ∂¯ t (ρu) + ∂ ∂¯ x ` ρu2 + p − Dρu ´ = 0 ∂ ∂¯ t „ ρ „ e + u2 2 «« + ∂ ∂¯ x „ u „ ρ „ e + u2 2 « + p « − Dρ „ e + u2 2 «« = 0 ∂ ∂¯ t (ρλ) + ∂ ∂¯ x (ρλu − Dρλ) = M ˙ ω

Shock Change Equation

At the shock, ρu = f(D) according to the jump conditions.

d dt(ρu)

˛ ˛

x=s

= dD

dt d dD (ρu)

˛ ˛

x=s

=

∂ ∂t(ρu)

˛ ˛

x=s + D ∂ ∂x(ρu)

˛ ˛

x=s .

Also the momentum equation gives

∂ ∂t(ρu) ˛ ˛ ˛ ˛

x=s

= − ∂ ∂x(ρu2 + p) ˛ ˛ ˛ ˛

x=s

.

Combining the two gives

dD dt = „ d dD (ρu) «−1 „ ∂ ∂x ` Dρu − ρu2 − p ´«˛ ˛ ˛ ˛ ˛

x=s

.

Numerics: WENO5 Mapped

x Stencil 2 Stencil 1

j j+1

Stencil 0

j+1/2

f

d f dx ˛ ˛ ˛ ˛

xj

= ˆ fj+1/2 − ˆ fj−1/2 ∆x

and

ˆ fj±1/2 =

2

X

k=0

ω(M)

k

ˆ f k

j±1/2,

where ω(M)

k

approximates the ideal weights to high order in smooth region of the flow.

• Conservative numerical scheme.
• Fifth order even near critical points of order one.
• Points near shock must still be considered.

Numerics: Boundary Points

x

Nx

ρ( ), D u(D), P(D)

x x

4 3

f ’(x) + O( ) f ’(x) + O( )

• Third and fourth order one-

sided derivatives at bound- ary points.

• Spatial order of the scheme

is at best third order globally.

• At the shock, conserved variables are functions of D

and the quiescent state.

• ∂

∂x at x = s for shock change equation calculated to

third order.

Numerics: Lax-Friedrichs Flux To properly discretize this system which involves both left and right traveling waves, use the Local Lax-Friedrichs flux:

f ± = f ± αu

where

α = |u − D| + c,

the maximum wave speed magnitude. The numerical flux ˆ

f is then found using WENO5M for both the f + and f − LLF: ˆ f = 1 2( ˆ f + + ˆ f −)

Numerics: Runge-Kutta Time Integration Integration in time using third-order TVD Runge-Kutta:

u∗ = un + ∆tL(un), u∗∗ = 3

4un + 1 4u∗ + 1 4∆tL(u∗),

un+1 = 1

3un + 2 3u∗∗ + 2 3∆tL(u∗∗),

where L = −

ˆ fj+1/2− ˆ fj+1/2 ∆x

.

ZND Structure Problem

D shocked state quiescent state reaction

0.2 0.4 0.6 0.8 1 1Ρ 10 20 30 40 50 60 70 p

• parameters:

E = 50, q = 50, γ = 1.2, Dcj = 6.80947

• quiescent state: ρ = 1, u = 0, p = 1, λ = 0
• consider f = (D/Dcj)2 = 1.8 and 1.6

9.13578 9.1358 9.13582 9.13584 9.13586 9.13588 9.1359 9.13592 9.13594 9.13596 20 40 60 80 100 D t Linear Decay for f=1.8 dx=0.1 dx=0.05

8 8.5 9 9.5 10 10.5 20 40 60 80 100 D t Non-Linear Saturation for f=1.6 dx=0.1 linear mode envelope

8.56 8.58 8.6 8.62 8.64 8.66 8.68 8.7 8.72 5 10 15 20 25 30 35 40 45 D t Linear Growth for f=1.6 dx=0.1 linear theory envelope

• 2.5
• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 7.5 8 8.5 9 9.5 10 10.5 D’(t) D Limit Cycle for f=1.6 dx=0.1

Conclusions

• Shock fitting allows for highly accurate solutions

containing a discontinuity.

• Validation of unstable mode for ZND problem from

linear theory.

• Explicit stability bounds and convergence study

presently lacking.

• 2-D extension needed for future work.