Effects of Diffusion on the Dynamics of Detonation Tariq D. Aslam - PowerPoint PPT Presentation

Effects of Diffusion on the Dynamics of Detonation Tariq D. Aslam Los Alamos National Laboratory; Los Alamos, NM Joseph M. Powers ∗ University of Notre Dame; Notre Dame, IN Gordon Research Conference - Energetic Materials Tilton, New Hampshire June 13-18, 2010 ∗ and contributions from many others

Motivation • Computational tools are critical in modeling of high speed reactive flow. • Steady wave calculations reveal sub-micron scale struc- tures in detonations with detailed kinetics (Powers and Paolucci, AIAA J. , 2005). • Small structures are continuum manifestation of molec- ular collisions. • We explore the transient behavior of detonations with fully resolved detailed kinetics.

Verification and Validation • verification: solving the equations right (math). • validation: solving the right equations (physics). • Main focus here on verification • Some limited validation possible, but detailed valida- tion awaits more robust measurement techniques. • Verification and validation always necessary but never sufficient: finite uncertainty must be tolerated.

Some Length Scales inherent in PBXs Micrograph of PBX 9501 (from C. Skidmore)

Some Length Scales Due to Diffusion Shock Rise in Aluminum (from V. Whitley) 10 ps rise time at 10 km/s yields scale of 10 − 7 m .

Modeling Issues for PBXs: • Inherently 3D, multi-component mixture, • Massive disparity in scales, • Many parameters are needed and many are unknown ∗ : – elastic constants – equation of states – thermal conductivities, viscosities for constituents – heat capacities – reaction rates – species diffusion ∗ see Menikoff and Sewell, CTM 2002

Before climbing Everest, we need to step back a bit... Let’s examine detonation dynamics of gases... 1. Inviscid, one-step Arrhenius chemistry 2. Inviscid, detailed chemistry 3. Diffusive, one-step Arrhenius chemistry 4. Diffusive, detailed chemistry

Model: Reactive Euler Equations • one-dimensional, • unsteady, • inviscid, • detailed mass action kinetics with Arrhenius tempera- ture dependency, • ideal mixture of calorically imperfect ideal gases

General Review of Pulsating Detonations • Erpenbeck, Phys. Fluids , 1962, • Fickett and Wood, Phys. Fluids , 1966, • Lee and Stewart, JFM , 1990, • Bourlioux, et al. , SIAM J. Appl. Math. , 1991, • He and Lee, Phys. Fluids , 1995, • Short, SIAM J. Appl. Math. , 1997, • Sharpe, Proc. R. Soc. , 1997.

Review of Recent Work of Special Relevance • Kasimov and Stewart, Phys. Fluids , 2004: published detailed discussion of limit cycle behavior with shock- fitting; error ∼ O (∆ x ) . • Ng, Higgins, Kiyanda, Radulescu, Lee, Bates, and Nikiforakis, CTM , in press, 2005: in addition, consid- ered transition to chaos; error ∼ O (∆ x ) . • Present study similar to above, but error ∼ O (∆ x 5 ) .

Model: Reactive Euler Equations • one-dimensional, • unsteady, • inviscid, • one step kinetics with finite activation energy, • calorically perfect ideal gases with identical molecular masses and specific heats.

Model: Reactive Euler Equations ∂ρ ∂t + ∂ ∂ξ ( ρu ) = 0 , ∂t ( ρu ) + ∂ ∂ ` ρu 2 + p ´ = 0 , ∂ξ ∂ „ „ e + 1 «« + ∂ „ „ e + 1 2 u 2 + p «« 2 u 2 = 0 , ρ ρu ∂t ∂ξ ρ „ « ∂t ( ρλ ) + ∂ ∂ − ρE ∂ξ ( ρuλ ) = αρ (1 − λ ) exp , p 1 p e = ρ − λq. γ − 1

Unsteady Shock Jump Equations ρ s ( D ( t ) − u s ) = ρ o ( D ( t ) − u o ) , � 1 � − 1 p s − p o = ( ρ o ( D ( t ) − u o )) 2 , ρ o ρ s � 1 e s − e o = 1 − 1 � 2( p s + p o ) , ρ o ρ s λ s = λ o .

Model Refinement • Transform to shock attached frame via � t x = ξ − D ( τ ) dτ, 0 • Use jump conditions to develop shock-change equa- tion for shock acceleration: � − 1 � ∂ � d ( ρ s u s ) � dD dt = − ∂x ( ρu ( u − D ) + p ) . dD

Numerical Method • point-wise method of lines, • uniform spatial grid, • fifth order spatial discretization (WENO5M) takes PDEs into ODEs in time only, • fifth order explicit Runge-Kutta temporal discretization to solve ODEs. • details in Henrick, Aslam, Powers, JCP , 2006.

Numerical Simulations • ρ o = 1 , p o = 1 , L 1 / 2 = 1 , q = 50 , γ = 1 . 2 , • Activation energy, E , a variable bifurcation parameter, 25 ≤ E ≤ 28 . 4 , √ � 61 • CJ velocity: D CJ = 11 + 5 ≈ 6 . 80947463 , • from 10 to 200 points in L 1 / 2 , • initial steady CJ state perturbed by truncation error, • integrated in time until limit cycle behavior realized.

Stable Case, E = 25 : Kasimov’s Shock-Fitting • N 1 / 2 = 100 , 200 , 6.84 N = 100 1/2 6.83 • minimum error in D : D N = 200 1/2 6.82 ∼ 9 . 40 × 10 − 3 , exact 6.81 • Error in D converges 0 50 100 150 200 250 300 t at O (∆ x 1 . 01 ) .

Stable Case, E = 25 : Improved Shock-Fitting • N 1 / 2 = 20 , 40 , 6.809485 • minimum error in D : 6.809480 ∼ 6 . 00 × 10 − 8 , for D exact 6.809475 N 1 / 2 = 40 . N = 20 1/2 6.809470 • Error in D converges 6.809465 0 50 100 150 200 250 300 at O (∆ x 5 . 01 ) . t

Linearly Unstable, Non-linearly Stable Case: E = 26 • One linearly unstable mode, stabilized by 7.4 non-linear effects, 7.0 D • Growth rate and fre- 6.6 quency match linear 6.2 0 100 200 300 400 500 theory to five decimal t places.

D, dD dt Phase Plane: E = 26 • Unstable spiral at early 0.4 time, stable period-1 stationary limit cycle 0.2 limit cycle at late time, dD dt 0.0 • Bifurcation point of E = 25 . 265 ± 0 . 005 -0.2 6.2 6.6 7.0 7.4 agrees with linear D stability theory.

Period Doubling: E = 27 . 35 • N 1 / 2 = 20 , 8 • Bifurcation to period- D 7 2 oscillation at E = 27 . 1875 ± 0 . 0025 . 6 0 100 200 300 400 t

D, dD dt Phase Plane: E = 27 . 35 stationary • Long time period-2 1 limit cycle limit cycle, dD dt 0 • Similar to independent results of Sharpe and -1 6 7 8 D Ng.

Transition to Chaos and Feigenbaum’s Number n →∞ δ n = E n − E n − 1 lim = 4 . 669201 . . . . E n +1 − E n E n +1 − E i n E n δ n 25 . 265 ± 0 . 005 0 - - 27 . 1875 ± 0 . 0025 1 . 9225 ± 0 . 0075 3 . 86 ± 0 . 05 1 2 27 . 6850 ± 0 . 001 0 . 4975 ± 0 . 0325 4 . 26 ± 0 . 08 3 27 . 8017 ± 0 . 0002 0 . 1167 ± 0 . 0012 4 . 66 ± 0 . 09 4 27 . 82675 ± 0 . 00005 0 . 02505 ± 0 . 00025 - . . . . . . . . . . . . ∞ 4 . 669201 . . .

Bifurcation Diagram 27.1875 < E < 27.6850 Stable Period 5 Stable Period 2 n 10 Period 2 1 9 E < 25.265, linearly stable 25.265 < E < 27.1875, D max 0 Period 2 (single period mode) 8 Stable 7 Period 6 Stable Period 3 25 26 27 28 E

D versus t for Increasing E 9 9 a b 8 8 D D 7 7 6 6 3700 3800 3900 4000 3700 3800 3900 4000 t t 9 c d 9 8 8 D D 7 7 6 6 3700 3800 3900 4000 3700 3800 3900 4000 t t 10 e f 9 9 8 8 D D 7 7 6 6 3700 3800 3900 4000 3700 3800 3900 4000 t t

Model: Reactive Euler PDEs with Detailed Kinetics ∂ρ ∂t + ∂ ∂x ( ρu ) = 0 , ∂t ( ρu ) + ∂ ∂ ρu 2 + p � � = 0 , ∂x e + u 2 e + u 2 � � �� � � �� ∂ + ∂ 2 + p ρ ρu = 0 , ∂t 2 ∂x ρ ∂t ( ρY i ) + ∂ ∂ ∂x ( ρuY i ) = M i ˙ ω i , N Y i � ρ ℜ T p = , M i i =1 e = e ( T, Y i ) , ω i ˙ = ω i ( T, Y i ) . ˙

Computational Methods • Steady wave structure – LSODE solver with IMSL DNEQNF for root finding – Ten second run time on single processor machine. – see Powers and Paolucci, AIAA J. , 2005. • Unsteady wave structure – Shock fitting coupled with a high order method for continuous regions – see Henrick, Aslam, Powers, J. Comp. Phys. , 2006, for full details on shock fitting

Ozone Reaction Kinetics a f β f E f j , a r j , β r j , E r Reaction j j j 6 . 76 × 10 6 1 . 01 × 10 12 O 3 + M ⇆ O 2 + O + M 2 . 50 1 . 18 × 10 2 3 . 50 0 . 00 4 . 58 × 10 6 2 . 51 × 10 11 O + O 3 ⇆ 2 O 2 2 . 50 1 . 18 × 10 6 4 . 15 × 10 12 2 . 50 5 . 71 × 10 6 4 . 91 × 10 12 O 2 + M ⇆ 2 O + M 2 . 50 2 . 47 × 10 2 3 . 50 0 . 00 see Margolis, J. Comp. Phys. , 1978, or Hirschfelder, et al. , J. Chem. Phys. , 1953.

Validation: Comparison with Observation • Streng, et al. , J. Chem. Phys. , 1958. • p o = 1 . 01325 × 10 6 dyne/cm 2 , T o = 298 . 15 K , Y O 3 = 1 , Y O 2 = 0 , Y O = 0 . Value Streng, et al. this study 1 . 863 × 10 5 cm/s 1 . 936555 × 10 5 cm/s D CJ T CJ 3340 K 3571 . 4 K 3 . 1188 × 10 7 dyne/cm 2 3 . 4111 × 10 7 dyne/cm 2 p CJ Slight overdrive to preclude interior sonic points.

Stable Strongly Overdriven Case: Length Scales D = 2 . 5 × 10 5 cm/s . −3 10 −4 10 length scale (cm) −5 10 −6 10 −7 10 −8 10 −10 −8 −6 −4 −2 0 10 10 10 10 10 10 x (cm)

Mean-Free-Path Estimate • The mixture mean-free-path scale is the cutoff mini- mum length scale associated with continuum theories. • A simple estimate for this scale is given by Vincenti and Kruger, ’65 : M 2 N πd 2 ρ ∼ 10 − 7 cm. √ ℓ mfp =

Stable Strongly Overdriven Case: Mass Fractions D = 2 . 5 × 10 5 cm/s . 10 1 O 2 10 0 O 10 −1 10 −2 Y i 10 −3 O 10 −4 3 −5 10 −9 −8 −7 −6 −5 −4 −3 −2 10 10 10 10 10 10 10 10 x (cm)