Resolution Matters: Issues in Computational Simulation of Detailed - PowerPoint PPT Presentation

Resolution Matters: Issues in Computational Simulation of Detailed Kinetics Gas Phase Combustion Joseph M. Powers (powers@nd.edu) Samuel Paolucci (paolucci@nd.edu) University of Notre Dame Notre Dame, Indiana 58 th Annual Meeting of the APS DFD Chicago, Illinois; 20 November 2005

Motivation • Detailed kinetics models are widely used in detonation simulations. • The finest length scale predicted by such models is usually not clarified and often not resolved. • Tuning computational results to match experiments without first harmonizing with underlying mathematics renders predictions unreliable. • See Powers and Paolucci, AIAA Journal , 2005.

Model: Steady 1D Reactive Euler Equations ρu = ρ o D, ρu 2 + p = ρ o D 2 + p o , e + u 2 ρ = e o + D 2 2 + p 2 + p o ρ o , N Y i � p = ρ ℜ T M i , i =1 � T � � N T − ℜ T � c pi ( ˆ T ) d ˆ h o e = Y i i,f + , M i T o i =1   � ρY k � ρY k � − Ej � ν ′ � ν ′′   J N N � − 1 dY i dx = M i kj kj � � �   ν ij α j T β j e ℜ T   K c ρ o D M k M k   j j =1   k =1 k =1 � �� � � �� � forward reverse

Eigenvalue Analysis of Local Length Scales Algebraic reduction yields d Y dx = f ( Y ) . Local behavior is modeled by d Y dx = J · ( Y − Y ∗ ) + b , Y ( x ∗ ) = Y ∗ . whose solution is Y ( x ) = Y ∗ + � � P · e Λ ( x − x ∗ ) · P − 1 − I · J − 1 · b . Here, Λ has eigenvalues λ i of Jacobian J in its diagonal. Length scales given by 1 ℓ i ( x ) = | λ i ( x ) | .

Computational Methods • A standard ODE solver (DLSODE) was used to inte- grate the equations. • Standard IMSL subroutines were used to evaluate the local Jacobians and eigenvalues at every step. • The Chemkin software package was used to evaluate kinetic rates and thermodynamic properties. • Computation time was typically one minute on a 1 GHz HP Linux machine.

Physical System • Hydrogen-air detonation: 2 H 2 + O 2 + 3 . 76 N 2 . • N = 9 molecular species, L = 3 atomic elements, J = 19 reversible reactions. • p o = 1 atm . • T o = 298 K . • Identical to system studied by both Shepherd (1986) and Mikolaitis (1987).

Mole Fractions versus Distance • significant evolution at 0 10 fine length scales x < N 2 H 2 H 2 H O O 2 2 OH −2 10 10 − 3 cm . O 2 O H −4 10 HO 2 H O −6 10 2 2 X i OH • results agree with −8 10 O −10 10 H O 2 2 those of Shepherd. −12 10 −14 10 −5 −4 −3 −2 −1 0 1 10 10 10 10 10 10 10 x (cm)

Eigenvalue Analysis: Length Scale Evolution • Finest length scale: 2 . 3 × 10 − 5 cm . 2 10 • Coarsest length scale 1 10 3 . 0 × 10 1 cm . 0 10 −1 10 (cm) i −2 10 • Finest length scale −3 10 similar to that −4 10 −5 10 necessary for −5 −4 −3 −2 −1 0 1 10 10 10 10 10 10 10 x (cm) numerical stability of ODE solver.

Numerical Stability • Discretizations finer than finest physical length −5 10 scale are numerically X H −6 stable. 10 -4 ∆ x = 2.38 x 10 cm (unstable) • Discretizations coarser -4 ∆ x = 2.00 x 10 cm (stable) −7 10 -5 ∆ x = 1.00 x 10 cm (stable) −4 −3 −2 than finest physical 10 10 10 x (cm) length scale are numerically unstable.

Examination of Recently Published Results ℓ ind ( cm ) ℓ f ( cm ) ∆ x ( cm ) Reference Under-resolution 2 × 10 − 1 2 × 10 − 4 4 × 10 − 3 2 × 10 1 Oran, et al. , 1998 2 × 10 − 2 5 × 10 − 5 3 × 10 − 3 6 × 10 1 Jameson, et al. , 1998 2 × 10 − 2 1 × 10 − 5 5 × 10 − 4 5 × 10 1 Hayashi, et al. , 2002 2 × 10 − 1 2 × 10 − 4 3 × 10 − 3 2 × 10 1 Hu, et al. , 2004 2 × 10 − 2 3 × 10 − 5 8 × 10 − 5 3 × 10 0 Powers, et al. , 2001 2 × 10 − 2 3 × 10 − 5 3 × 10 − 2 1 × 10 3 Osher, et al. , 1997 5 × 10 − 3 8 × 10 − 6 1 × 10 − 2 1 × 10 3 Merkle, et al. , 2002 1 × 10 − 1 2 × 10 − 4 1 × 10 0 5 × 10 3 Sislian, et al. , 1998 2 × 10 − 2 6 × 10 − 7 6 × 10 − 2 1 × 10 5 Jeung, et al. , 1998 All are under-resolved, some severely.

Conclusions • Detonation calculations are often under-resolved, by as much as five orders of magnitude. • Equilibrium properties are insensitive to resolution, while transient phenomena can be sensitive. • Sensitivity of results to resolution is not known a priori . • Numerical viscosity stabilizes instabilities. • For a repeatable scientific calculation of detonation, the finest physical scales must be resolved.

Moral You either do detailed kinetics with the proper resolution, or you are fooling yourself and others, in which case you should stick with reduced kinetics!

Detailed Kinetics Model j Aj βj Ej Reaction 1 . 70 × 1013 H 2 + O 2 ⇀ ↽ OH + OH 0 . 00 47780 1 1 . 17 × 109 OH + H 2 ⇀ ↽ H 2 O + H 1 . 30 3626 2 5 . 13 × 1016 H + O 2 ⇀ ↽ OH + O − 0 . 82 16507 3 1 . 80 × 1010 O + H 2 ⇀ ↽ OH + H 1 . 00 8826 4 2 . 10 × 1018 H + O 2 + M ⇀ ↽ HO 2 + M − 1 . 00 5 0 6 . 70 × 1019 H + O 2 + O 2 ⇀ ↽ HO 2 + O 2 − 1 . 42 0 6 6 . 70 × 1019 H + O 2 + N 2 ⇀ ↽ HO 2 + N 2 − 1 . 42 7 0 5 . 00 × 1013 OH + HO 2 ⇀ ↽ H 2 O + O 2 0 . 00 1000 8 2 . 50 × 1014 H + HO 2 ⇀ ↽ OH + OH 0 . 00 1900 9 4 . 80 × 1013 O + HO 2 ⇀ ↽ O 2 + OH 0 . 00 1000 10 6 . 00 × 108 OH + OH ⇀ ↽ O + H 2 O 1 . 30 0 11 2 . 23 × 1012 H 2 + M ⇀ ↽ H + H + M 0 . 50 92600 12 1 . 85 × 1011 O 2 + M ⇀ ↽ O + O + M 0 . 50 95560 13 7 . 50 × 1023 H + OH + M ⇀ ↽ H 2 O + M − 2 . 60 14 0 2 . 50 × 1013 H + HO 2 ⇀ ↽ H 2 + O 2 0 . 00 700 15 2 . 00 × 1012 HO 2 + HO 2 ⇀ ↽ H 2 O 2 + O 2 0 . 00 16 0 1 . 30 × 1017 H 2 O 2 + M ⇀ ↽ OH + OH + M 0 . 00 45500 17 1 . 60 × 1012 H 2 O 2 + H ⇀ ↽ HO 2 + H 2 0 . 00 3800 18 1 . 00 × 1013 H 2 O 2 + OH ⇀ ↽ H 2 O + HO 2 0 . 00 1800 19

Temperature Profile • Temperature flat in the post-shock induction 3500 zone 0 < x < 3000 2 . 6 × 10 − 2 cm . T (K) 2500 • Thermal explosion 2000 followed by relaxation 1500 −5 −4 −3 −2 −1 0 1 10 10 10 10 10 10 10 to equilibrium at x (cm) x ∼ 10 0 cm .

Verification: Comparison with Mikolaitis • Lagrangian calculation −1 10 allows direct −3 10 −5 comparison with 10 X i HO 2 −7 10 Mikolaitis’ results. H −9 OH 10 O H O • agreement very good. 2 2 −11 10 H O 2 −13 10 −12 −11 −10 −9 −8 −7 −6 10 10 10 10 10 10 10 t (s)

Grid Convergence • Finest length scale must be resolved to converge at proper 0 10 −2 10 order. 1.006 First Order −4 10 Explicit Euler ε OH 1 • Results are −6 10 2.008 −8 converging at proper 10 Second Order Runge-Kutta 1 −10 10 order for first and −10 −8 −6 −4 10 10 10 10 ∆ x (cm) second order discretizations.