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 43 rd 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 y β Y u s s 1 i n o β c • Straight shock. 1 u u 1 p = p 1 reaction zone 1 ~ = ρ = ρ λ k • Curved wedge. streamline c 1 o h T = T s 1 t λ = 0 h g i a r t s • Orthogonal coordinate curved wedge β system aligned with x shock. X

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 , ∂ + ∂ ρU 2 + p ` ´ ∂Y ( ρUV ) = 0 , ∂X ∂X ( ρUV ) + ∂ ∂ ρV 2 + p ` ´ = 0 , ∂Y „ „ «« ∂ e + 1 2( U 2 + V 2 ) + p ρU ∂X ρ „ „ «« + ∂ e + 1 2( U 2 + V 2 ) + p ρV = 0 , ∂Y ρ ∂X ( ρUλ ) + ∂ ∂ ∂Y ( ρV λ ) = αρ (1 − λ ) H ( T − T i ) , 1 p = e ρ − λq, γ − 1 p = ρRT.

Model Reductions: PDEs → ODEs Assume no Y variation, so d dX ( ρU ) = 0 , d ρU 2 + p � � = 0 , dX d dX ( ρUV ) = 0 , d � � e + 1 2( U 2 + V 2 ) + p �� ρU = 0 , dX ρ d dX ( ρUλ ) = αρ (1 − λ ) H ( T − T i ) .

Model Reductions: ODEs → DAEs ρU = ρ 1 u 1 sin β, ρU 2 + p 1 sin 2 β + p 1 , ρ 1 u 2 = V = u 1 cos β, γ p ρ − λq + 1 γ p 1 + 1 U 2 + u 2 1 cos 2 β 2 u 2 � � = 1 , γ − 1 2 γ − 1 ρ 1 dλ α 1 − λ = H ( T − T i ) . dX U ZND reaction zone structure ODE supplemented with extended Rankine-Hugoniot algebraic conditions.

Model Reductions: Inversion of Algebraic Relations with M 1 ≡ M 1 sin β , ρ 1 ( γ + 1) M 2 1 ρ ( λ ) = ” , r` 1 + γ M 2 ´ 2 − ( γ + 1) M 2 “ 2 + γ − 1 2 λq 1 + γ M 2 RT 1 + ( γ − 1) M 2 1 ± 1 1 1 γ ρ 1 u 1 sin β U ( λ ) = , ρ ( λ ) „ 1 1 « 1 sin 2 β p 1 + ρ 2 1 u 2 p ( λ ) = − , ρ ( λ ) ρ 1 „ 1 1 sin 2 β ρ ( λ ) R + ρ 2 1 u 2 1 « p 1 T ( λ ) = − , ρ ( λ ) R ρ ( λ ) ρ 1 γRT 1 ( M 2 1 − 1) 2 CJ limitation . q ≤ , 2( γ 2 − 1) M 2 1 + shocked; − unshocked. Take the shocked branch.

Reaction Zone Structure Solution dλ α dX = ρ 1 u 1 sin β ρ ( λ )(1 − λ ) , λ (0) = 0 , p 1 − a 4 0 0 a 3 1 1 r 0 „ « „ « 1 1 − a 4 λ r 1 1 − 1 + ! a 2 1 1 − a 4 1 − a 4 B B C C “p ” B C B B C C X ( λ ) = a 1 2 a 3 1 − a 4 λ − 1 + ln , B C B B C C B „ r « „ « C 1 − a 4 λ r B B 1 − λ C C 1 @ A 1 + 1 − @ @ A A 1 − a 4 1 − a 4 √ γRT 1 1 a 1 = , ( γ + 1) M 1 α 1 + γ M 2 a 2 = 1 , M 2 a 3 = 1 − 1 , γ 2 − 1 M 2 q 1 a 4 = 2 . 1 − 1) 2 γ RT 1 ( M 2

Parametric Values Independent Parameter Units Value R J/kg/K 287 α 1 /s 1000 β rad π/ 4 γ 6 / 5 - T 1 K 300 M 1 3 - kg/m 3 ρ 1 1 q J/kg 300000 T i K 131300 / 363

Reaction Zone Structure Normal to Shock 3 λ ρ (kg/ m ) 1 4 0.8 a) 3 b) 0.6 2 0.4 1 0.2 X (m) 4 X (m) 0.0 0 -1 0 1 2 3 4 -1 0 1 2 3 T(K) MX 600 3 500 c) 2 d) 400 1 300 X (m) X (m) 200 0 2 -1 0 1 2 3 4 -1 0 1 3 4

Exact Solution: Streamlines y (m) straight oblique shock, β = π / 4 4 • Curved streamlines identical to wedge 3 contour. 2 λ ~ 0.9 • Streamline curvature 1 approaches zero as curved wedge 0 x (m) reaction completes. 0 1 2 4 3

High Mach Number Limit Solution „ 1 „ „ « «« ρ ( λ ) γ + 1 1 γ + 1 2 γ γ 2 − 1 + λq = 1 − + O , M 2 M 4 ρ 1 γ − 1 γ RT 1 1 1 „ 1 „ „ « «« U ( λ ) γ − 1 1 γ + 1 2 γ γ 2 − 1 + λq = 1 + + O , M 2 M 4 u 1 sin β γ + 1 γ RT 1 1 1 „ 1 γ 2 − 1 „ „ « «« p ( λ ) 2 γ 1 γ + 1 + λq 1 M 2 = 1 − + O . 1 γ + 1 M 2 2 γ M 4 p 1 RT 1 1 1 „ « ( γ + 1) α λ ( X ) = 1 − exp − ( γ − 1) u 1 sin β X , X r ∼ γ − 1 u 1 sin β Reaction zone thickness . γ + 1 α

Low Heat Release Limit Effects of heat release are better represented following a detailed asymptotic analysis, which yields a 5 ln(1 − λ ) − a 3 a 4 “ ” X ( λ ) = a 1 λ . 2 Invert to form „ X » a 3 a 4 «– λ ( X ) = 1 − 2 a 5 a 1 a 5 + a 3 a 4 2 a 5 exp a 3 a 4 W 0 . 2 a 5 „ « X r ∼ γ − 1 u 1 sin β 2 γ + 1 q 1 + + . γ + 1 ( γ − 1) M 2 γ ( M 2 1 − 1) α RT 1 1 Lambert W 0 function utilized: W 0 ( we w ) = w.

Exact versus Asymptotic Solutions • High Mach number λ 1.0 limit solution agrees exact solution 0.8 low heat release limit high Mach number limit poorly. 0.6 0.4 • Low heat release limit 0.2 X (m) 0.0 0 1 2 3 4 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 y (m) 2.0 • 256 × 256 uniform ρ = 2.9 kg/m 3 exact ρ = 2.6 kg/m 3 1.5 numerical ρ = 2.3 kg/m 3 numerical grid. 1.0 ρ = 2.0 kg/m 3 • good agreement in pic- 0.5 ture norm. 0.0 x (m) 0.0 0.5 1.0 1.5 • numerical solution sta- ble.

Iterative Convergence to Steady State: Various Grids • Coarse grids relax L (kg/m) 1 quickly; fine grids relax 10 64 x 64 slowly. 128 x 128 256 x 256 1 512 x 512 1024 x 1024 • All grids iteratively 0.1 converge to steady 0.01 state. t (s) 0 0.002 0.004 0.006 0.008 0.01 • Iterative convergence is distinct from grid convergence.

Grid Convergence L (kg/m) 1 • Convergence 1 rate: ∆ x 0 . 779 � � O . 0.1 • Both shock capturing and embedded bound- 0.01 ∆ x (m) 0.001 0.01 0.1 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 order convergence rates even for high order discretiza- tions. • Common practice of claiming high order convergence rates without verification should be stopped.