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 ) + ∂ ∂ ρu 2 + p � � = 0 ∂x e + u 2 e + u 2 ∂ � � �� + ∂ � � � � ρ u ( ρ + p ) = 0 ∂t 2 ∂x 2 � − Eρ � ∂t ( ρλ ) + ∂ ∂ ∂x ( ρλu ) = kρ (1 − λ ) exp p 1 p = ρRT e = γ − 1 p/ρ − λq

Jump Conditions Governing PDE’s do not hold over a discontinuity - σ (1) must appeal to integral formulation. Σ γ (1) Z Z Z σ (2) d V = + Fd Bd V T i n i dS − − dt M − V M − V MS γ (2) | {z } | {z } ν ι surface forces volume forces V (1) � F ( D i − v i ) + T i � = 0 V (2) � ρ ( D − u ) � = 0 � ρu ( D − u ) − p � = 0 „ « e + u 2 � ρ ( D − u ) − up � = 0 2 � ρλ ( D − u ) � = 0 This gives ρ ( D ) , u ( D ) , p ( D ) , λ ( D ) given the quiescent state.

Wave Frame Transformation = ∂ ¯ ∂t + ∂ ¯ ∂ ¯ ∂f f ∂ ¯ x f t 9 ∂ ¯ ∂t ∂ ¯ x t ∂t x ( x, t ) = x − s ( t ) ¯ −∞ < x < + ∞ = = ∂ ¯ s ( t )) + ∂ ¯ f f x ( − ˙ ∂ ¯ ¯ ∂ ¯ t t ( x, t ) = t all x and t, ; = − D ∂ ¯ x + ∂ ¯ f f ∂ ¯ ∂ ¯ t = ∂ ¯ ∂x + ∂ ¯ ∂ ¯ ∂f f f ∂ ¯ x t ∂ ¯ ∂x ∂ ¯ x t ∂x = ∂ ¯ f ∂ ¯ x ∂ρ t + ∂ x ( ρu − Dρ ) = 0 ∂ ¯ ∂ ¯ ∂ t ( ρu ) + ∂ ρu 2 + p − Dρu ` ´ = 0 ∂ ¯ ∂ ¯ x „ „ «« „ „ „ « « „ «« e + u 2 e + u 2 e + u 2 ∂ + ∂ + p = 0 ρ u ρ − Dρ ∂ ¯ t 2 ∂ ¯ x 2 2 ∂ t ( ρλ ) + ∂ x ( ρλu − Dρλ ) = M ˙ ω ∂ ¯ ∂ ¯

Shock Change Equation At the shock, ρu = f ( D ) according to the jump conditions. ˛ ˛ d = dD d dt ( ρu ) dD ( ρu ) ˛ ˛ dt x = s x = s ˛ ˛ ∂ ∂ = ∂t ( ρu ) x = s + D ∂x ( ρu ) x = s . ˛ ˛ Also the momentum equation gives ˛ ˛ ∂ = − ∂ ˛ ∂x ( ρu 2 + p ) ˛ ∂t ( ρu ) . ˛ ˛ ˛ ˛ x = s x = s Combining the two gives „ d « − 1 „ ∂ ´«˛ ˛ dD ` Dρu − ρu 2 − p ˛ dt = dD ( ρu ) . ˛ ∂x ˛ x = s

Numerics: WENO5 Mapped f ˛ f j +1 / 2 − ˆ ˆ j+1/2 f j − 1 / 2 d f ˛ = ˛ and x dx ∆ x ˛ x j j j+1 2 X ˆ ω ( M ) f k ˆ Stencil 0 f j ± 1 / 2 = j ± 1 / 2 , k k =0 Stencil 1 where ω ( M ) approximates the ideal weights to k Stencil 2 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 • Third and fourth order one- D u(D), P(D) ρ( ), x sided derivatives at bound- ary points. 0 Nx • Spatial order of the scheme f ’(x) + O( ) x 4 f ’(x) + O( ) x 3 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 α = | u − D | + c, where the maximum wave speed magnitude. The numerical flux ˆ f is then found using WENO5M for both the f + and f − LLF: f = 1 f + + ˆ ˆ 2( ˆ f − )

Numerics: Runge-Kutta Time Integration Integration in time using third-order TVD Runge-Kutta: u ∗ = u n + ∆ t L ( u n ) , u ∗∗ = 3 4 u n + 1 4 u ∗ + 1 4 ∆ t L ( u ∗ ) , u n +1 = 1 3 u n + 2 3 u ∗∗ + 2 3 ∆ t L ( u ∗∗ ) , f j +1 / 2 − ˆ ˆ f j +1 / 2 where L = − . ∆ x

ZND Structure Problem p shocked state 70 60 quiescent state 50 ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ �✁�✁�✁�✁�✁�✁�✁�✁� D 40 30 20 10 reaction 1 1 � Ρ 0.2 0.4 0.6 0.8 • parameters: E = 50 , q = 50 , γ = 1 . 2 , D cj = 6 . 80947 • quiescent state: ρ = 1 , u = 0 , p = 1 , λ = 0 • consider f = ( D/D cj ) 2 = 1 . 8 and 1 . 6

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

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

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

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

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.