Exact Solutions for Two-Dimensional Reactive Flow for Verification - - PowerPoint PPT Presentation

Exact Solutions for Two-Dimensional Reactive Flow for Verification of Numerical Algorithms Joseph M. Powers (powers@nd.edu) University of Notre Dame; Notre Dame, IN Tariq D. Aslam (aslam@lanl.gov) Los Alamos National Laboratory; Los Alamos, NM 43rd AIAA Aerospace Sciences Meeting and Exhibit Reno, Nevada 10-13 January 2005

Motivation

• Computational tools are critical in design of aerospace

vehicles which employ high speed reactive flow.

• Comparing computational predictions with those of

exact solutions in grid resolution studies is a robust verification.

• We develop a new exact solution and employ it to verify

a modern shock-capturing reactive flow algorithm for flows with an immersed boundary.

Verification and Validation

• verification: solving the equations right.
• validation: solving the right equations.
• Focus here is exclusively on verification.
• Limiting assumptions necessary for exact solution pre-

clude meaningful validation exercise.

• Verification and validation always necessary but never

sufficient: finite uncertainty must be tolerated.

Partial Review of Oblique Detonations

• Samaras, Can. J. Research , 1948.
• Gross, AIAA J., 1963.
• Lee, AIAA J., 1966.
• Pratt, J. Propul. Power, 1991.
• Powers, et al., AIAA J., Phys. Fluids, Shock Waves,

1992-96.

Oblique Detonation Schematic

x y

X Y

β s t r a i g h t s h

• c

k curved wedge

u1 u c

• s

β

1

u s i n β

1

p = p ρ = ρ T = T λ = 0

1 1 1 λ = 1

streamline reaction zone ~

• Straight shock.
• Curved wedge.
• Orthogonal coordinate

system aligned with shock.

Model: Reactive Euler Equations

• two-dimensional,
• steady,
• inviscid,
• irrotational,
• one step kinetics with zero activation energy,
• calorically perfect ideal gases with identical molecular

masses and specific heats.

Model: Reactive Euler PDEs

∂ ∂X (ρU) + ∂ ∂Y (ρV ) = 0, ∂ ∂X ` ρU 2 + p ´ + ∂ ∂Y (ρUV ) = 0, ∂ ∂X (ρUV ) + ∂ ∂Y ` ρV 2 + p ´ = 0, ∂ ∂X „ ρU „ e + 1 2(U 2 + V 2) + p ρ «« + ∂ ∂Y „ ρV „ e + 1 2(U 2 + V 2) + p ρ «« = 0, ∂ ∂X (ρUλ) + ∂ ∂Y (ρV λ) = αρ(1 − λ)H(T − Ti), e = 1 γ − 1 p ρ − λq, p = ρRT.

Model Reductions: PDEs→ODEs Assume no Y variation, so

d dX (ρU) = 0, d dX

• ρU 2 + p
• =

0, d dX (ρUV ) = 0, d dX

• ρU
• e + 1

2(U 2 + V 2) + p ρ

• =

0, d dX (ρUλ) = αρ(1 − λ)H(T − Ti).

Model Reductions: ODEs→DAEs

ρU = ρ1u1 sin β, ρU 2 + p = ρ1u2

1 sin2 β + p1,

V = u1 cos β, γ γ − 1 p ρ − λq + 1 2

• U 2 + u2

1 cos2 β

• =

γ γ − 1 p1 ρ1 + 1 2u2

1,

dλ dX = α1 − λ U H(T − Ti).

ZND reaction zone structure ODE supplemented with extended Rankine-Hugoniot algebraic conditions.

Model Reductions: Inversion of Algebraic Relations

with M1 ≡ M1 sin β, ρ(λ) = ρ1(γ + 1)M2

1

1 + γM2

1 ±

r`1 + γM2

1

´2 − (γ + 1)M2

1

“ 2 + γ−1

γ 2λq RT1 + (γ − 1)M2 1

”, U(λ) = ρ1u1 sin β ρ(λ) , p(λ) = p1 + ρ2

1u2 1 sin2 β

„ 1 ρ1 − 1 ρ(λ) « , T(λ) = p1 ρ(λ)R + ρ2

1u2 1 sin2 β

ρ(λ)R „ 1 ρ1 − 1 ρ(λ) « , q ≤ γRT1(M2

1 − 1)2

2(γ2 − 1)M2

1

,

CJ limitation.

+ shocked; − unshocked. Take the shocked branch.

Reaction Zone Structure Solution

dλ dX = α ρ1u1 sin β ρ(λ)(1 − λ), λ(0) = 0,

X(λ) = a1 B B B B @ 2a3 “p 1 − a4λ − 1 ” + ln B B B B @ 1 1 − λ !a2 B B B @ „ 1 − r 1−a4λ 1−a4 « „ 1 + r 1 1−a4 « „ 1 + r 1−a4λ 1−a4 « „ 1 − r 1 1−a4 « 1 C C C A a3 p1−a4 1 C C C C A 1 C C C C A ,

a1 = 1 (γ + 1)M1 √γRT1 α , a2 = 1 + γM2

1,

a3 = M2

1 − 1,

a4 = 2 M2

1

(M2

1 − 1)2

γ2 − 1 γ q RT1 .

Parametric Values

Independent Parameter Units Value

R J/kg/K 287 α 1/s 1000 β rad π/4 γ

• 6/5

T1 K 300 M1

• 3

ρ1 kg/m3 1 q J/kg 300000 Ti K 131300/363

Reaction Zone Structure Normal to Shock

• 1 0

1 2 3 4 X (m) 1 2 3 4 ρ (kg/ m )

3

• 1 0

1 2 3 4 X (m) 300 400 500 600 T(K)

• 1 0

1 2 3 4 X (m) 1 2 3 MX

• 1 0

1 2 3 4 X (m) 0.2 0.4 0.6 0.8 1 λ

a) b) c) d)

200 0.0

Exact Solution: Streamlines

2 4 1 2 3 4 1 3 x (m) y (m) λ ~ 0.9 straight oblique shock, β = π / 4 curved wedge

• Curved

streamlines identical to wedge contour.

• Streamline

curvature approaches zero as reaction completes.

High Mach Number Limit Solution

ρ(λ) ρ1 = γ + 1 γ − 1 „ 1 − 1 M2

1

γ + 1 γ „ 2γ γ2 − 1 + λq RT1 « + O „ 1 M4

1

«« , U(λ) u1 sin β = γ − 1 γ + 1 „ 1 + 1 M2

1

γ + 1 γ „ 2γ γ2 − 1 + λq RT1 « + O „ 1 M4

1

«« , p(λ) p1 = M2

1

2γ γ + 1 „ 1 − 1 M2

1

γ2 − 1 2γ „ 1 γ + 1 + λq RT1 « + O „ 1 M4

1

«« . λ(X) = 1 − exp „ − (γ + 1)α (γ − 1)u1 sin β X « , Xr ∼ γ − 1 γ + 1 u1 sin β α .

Reaction zone thickness

Low Heat Release Limit Effects of heat release are better represented following a detailed asymptotic analysis, which yields

X(λ) = a1 “ a5 ln(1 − λ) − a3a4 2 λ ” .

Invert to form

λ(X) = 1 − 2a5 a3a4 W0 »a3a4 2a5 exp „ X a1a5 + a3a4 2a5 «– . Xr ∼ γ − 1 γ + 1 u1 sin β α „ 1 + 2 (γ − 1)M2

1

+ γ + 1 γ(M2

1 − 1)

q RT1 « .

Lambert W0 function utilized:

W0(wew) = w.

Exact versus Asymptotic Solutions

1 2 3 4 0.2 0.4 0.6 0.8 1.0 0.0

X (m) λ

exact solution low heat release limit high Mach number limit

• High

Mach number limit solution agrees poorly.

• Low heat release limit

solution agrees well.

Verification of Modern Shock Capturing Algorithm

• Algorithm of Xu, Aslam, and Stewart, 1997, CTM.
• Uniform Cartesian grid.
• Embedded internal boundary with level set represen-

tation.

• Nominally fifth order weighted essentially non-oscillatory

(WENO) discretization.

• Non-decomposition based Lax-Friedrichs solver.
• Third order Runge-Kutta time integration.

Exact versus Numerical Solutions

0.0 0.5 1.0 1.5 0.0 0.5 1.0 1.5 2.0

x (m) y (m) exact numerical

ρ = 2.0 kg/m3 ρ = 2.9 kg/m3 ρ = 2.6 kg/m3 ρ = 2.3 kg/m3

• 256 × 256 uniform

numerical grid.

• good agreement in pic-

ture norm.

• numerical solution sta-

ble.

Iterative Convergence to Steady State: Various Grids

0.01 0.1 1 10 0.002 0.004 0.006 0.008 0.01 64 x 64 128 x 128 256 x 256 512 x 512 1024 x 1024 t (s) L (kg/m)

1

• Coarse

grids relax quickly; fine grids relax slowly.

• All

grids iteratively converge to steady state.

• Iterative

convergence is distinct from grid convergence.

Grid Convergence

0.01 0.1 1 0.001 0.01 0.1 ∆x (m) L (kg/m)

1

• Convergence

rate:

O

• ∆x0.779

.

• Both shock capturing

and embedded bound- ary induce the low con- vergence rate.

Conclusions

• New exact solution for two-dimensional steady detona-

tion found.

• Excellent verification tool for computational methods.
• Numerical solutions are stable.
• Shock capturing and embedded boundary induce low
• rder convergence rates even for high order discretiza-

tions.

• Common practice of claiming high order convergence

rates without verification should be stopped.