Numerical Approach on Hydrogen Numerical Approach on Hydrogen Detonation: Fundamentals and Detonation: Fundamentals and Applications Applications -Part 1 Part 1- - - 2007.08.02 2007.08.02 Nobuyuki TSUBOI Nobuyuki TSUBOI ISAS/JAXA, Japan ISAS/JAXA, Japan 1 1
Overview Overview 1.Motivations 2. Numerical simulation for compressible high speed flow (1) Scalar equations (2) System equations (3) System equations with chemical reactions 3.Summary 2 2
Motivations Motivations 1.Hydrogen/air mixture: detonable gas 2.Detonation: shock induced combustion -Pressure behind detonation increases about 10 times ambient pressure 3.Closed environment such as a tunnel causes serious accident. 3 3
Numerical Simulations Numerical Simulations 1.Detonation as well as thermal and gas dynamic phenomena of airplane and aerospace plane are elucidated numerically. 2. Main issues (1) Robust numerical scheme for high pressure and temperature (2) Stiff problem: chemical characteristic time is much faster than the fluid characteristic time (3) Chemical reaction model for high pressure combustion (4) Unsteadiness, turbulence, and three- dimensional problems 4 4
Flow Chart of Numerical Simulations Flow Chart of Numerical Simulations Flowfield to Simulate ? ・・ Compressible? Incompressible? Physical Model ? Chemical reaction? Turbulence? ・・ Euler equations ? Mathematical Model ? Navier-Stokes equations ? ・・ Finite difference ? Finite Volume ? Discretization ? Finite element ? Computational Grid ・・ Structure grid ? Unstructure grid ? Cartesian grid ? Program ・・ Fortran? C? ・・ Supercomputer? Workstation? Submit Job ・・ Workstation? Visualization Personal computer? 5 5
Numerical Model Numerical Model Level Mathematical Model Target Level Mathematical Model Target Panel Method w/o separation Conceptual Conceptual Panel Method w/o separation design design Potential Eq Eq. . less turbulence Potential less turbulence Conceptual Euler Eq Eq. . Shock Conceptual Euler Shock design design Academic RANS(Reynolds averaged averaged Boundary layer, Academic RANS(Reynolds Boundary layer, research, research, Navier- - Stokes Stokes Eq Eq.) .) separation Navier separation Detailed Detailed design design Academic LES(Large Eddy Simulation) Eddy Simulation) Solve scale Academic LES(Large Solve scale research research vortices smaller vortices smaller than grid than grid DNS(Direct Numerical Numerical Solve all vortices Academic Academic DNS(Direct Solve all vortices research research Simulation) Simulation) 6 6
Computational Grids Computational Grids Structured grid ? Unstructured grid ? Cartesian grid ? How to make grids ? ・・ Commercial software ? ・・ Make by oneself ? 7 7
Numerical analysis for fluid: finite different method Numerical analysis for fluid: finite different method Flowchart of program is defined after discretization and grid are decided Start Initialization Read grid, parameter, and set initial condition Boundary condition CFL number Define time step Convective flux, viscous term Upwind method, central difference, No Calculate numerical flux Higher-order Calculate time integration Explicit or implicit time integration Convergence ? Output data Output simulation results End 8 8
Scalar Equation Scalar Equation 9 9
Finite Difference for Scalar Equation Finite Difference for Scalar Equation Wave equation ・・ wave propagates with speed, c ∂ ∂ u u + = c 0 ∂ ∂ t x u c u 1 u 2 x 10 10
∂ ∂ u u + = c 0 ∂ ∂ t x FTCS (Forward in Time and Central Difference in Space) ・・ Explicit+central difference:unstable Δ ⎛ ⎞ c t + = − − ⎜ ⎟ n 1 n n n Δ ( ) ・・ CFL number: u u u u ⎛ ⎞ t + − Δ ν = j j j 1 j 1 ⎝ ⎠ ⎜ ⎟ c 2 x parameter to govern Δ ⎝ ⎠ x numerical stability Lax method Δ ( ) ⎛ ⎞ 1 c t More stable than FTCS + = + − − ⎜ ⎟ n 1 n n n n u u u ( u u ) − + + − Δ Large numerical dissipation j j 1 j 1 j 1 j 1 ⎝ ⎠ 2 2 x 2nd order in space and time Lax-Wendroff method Δ Δ 2 ⎛ ⎞ ⎛ ⎞ 2 ( ) ( ) c t c t + = − − + − + n 1 n n n n n n ⎜ ⎟ ⎜ ⎟ u u u u u 2 u u + − + − Δ Δ j j j 1 j 1 j 1 j j 1 ⎝ ⎠ ⎝ ⎠ 2 x 2 x First-order upwind method (Upwind difference) Δ ⎛ ⎞ t + = − − ⎜ ⎟ ・・ Find direction which wave propagates n 1 n n n u u c ( u u ) − Δ j j j j 1 ⎝ ⎠ x and use upwind data only 11 11
Pseudo-finite volume method ・ Conservation form ・・ Non-linear equation with discontinuous wave can be calculate and discreted ・・ Compressible equations is based on this concept ∂ ∂ u f + = = 0, f cu Flux ∂ ∂ t x Δ ⎛ ⎞ ( ) t ~ % % + = − − n n 1 n n n ⎜ ⎟ f u u f f : Numerical flux + + − Δ j 1 / 2 j j j 1/2 j 1/2 ⎝ ⎠ x − + 1 j 1 j 2 2 ~ % n n f − f + j 1 / 2 j 1/ 2 − + j j 1 j 1 12 12
Δ ⎛ ⎞ ( ) t % % + = − − 1 n n n n ⎜ ⎟ u u f f + − Δ j j j 1/2 j 1/2 ⎝ ⎠ x For example, FTCS is applied, Δ ⎛ ⎞ ( ) c t + = − − n 1 n n n ⎜ ⎟ u u u u + − Δ j j j 1 j 1 ⎝ ⎠ 2 x Δ { } ⎛ ⎞ ( ) ( ) 1 t = − + − + n n n n n ⎜ ⎟ u cu cu cu cu + − Δ j j 1 j j j 1 ⎝ ⎠ 2 x Δ ⎛ ⎞ ( ) t % % = − − n n n ⎜ ⎟ u f f + − Δ j ⎝ ⎠ j 1/2 j 1/ 2 x Therefore ( ) ( ) 1 c % = + = + n n n n n f f f u u + + + j 1/2 j 1 j j 1 j 2 2 13 13
There exist many numerical schemes, however, these schemes are designed ~ n how to estimate numerical flux f + j 1 / 2 FTCS n + f j = 1 n ) = c ˜ + u j n n n ) f 2 ( f j + 1 2 ( u j + 1 j + 1 / 2 Lax Δ ⎛ ⎞ ~ c 1 x = + − − ⎜ ⎟ n n n n n f ( u u ) ( u u ) + + + Δ j 1 / 2 j 1 j j 1 j ⎝ ⎠ 2 2 t Lax-Wendroff ~ c c 1 = + ν + − ν = + − ν − n f ( 1 ) u ( 1 ) u c [ u ( 1 )( u u )] + + + 1 / 2 1 1 j j j j j j 2 2 2 Δ ⎛ ⎞ 1st order upwind difference t ν = ⎜ ⎟ c Δ n = f j ⎝ ⎠ = c n ) − c x ˜ + u j − u j n ) = cu j n n n n f 2 ( u j + 1 2 ( u j + 1 j + 1 / 2 14 14
1st order upwind difference n = f j = c n ) − c ˜ + u j − u j n ) = cu j n n n n 2 ( u j + 1 2 ( u j + 1 f j + 1 / 2 When sign of c is unclear, summarized as follows: − + c c c c ~ = + n f u u + + 1 / 2 1 j j j ⎧ ≤ 2 2 ⎪ n f , 0 c ~ [ ] = j n ⎨ 1 f = + − − ( cu cu ) c ( u u ) + ≤ j 1 / 2 ⎪ + + n j 1 j j 1 j f , c 0 ⎩ 2 + j 1 [ ] 1 = + − − ( ) ( ) f f c u u + + j 1 j j 1 j 2 15 15
TVD(Total Variation Variation Dimishing Dimishing) Method ) Method TVD(Total 1st order upwind ˜ = cu j n n f j + 1 / 2 difference ~ c c Lax-Wendroff = + ν + − ν n f ( 1 ) u ( 1 ) u + + j 1 / 2 j j 1 2 2 1 = + − ν − c [ u ( 1 )( u u )] + j j 1 j 2 L-W scheme = 1st order upwind + modified flux -> 2nd order 16 16
~ 1 = + − ν − n [ ( 1 )( )] f c u u u + + j 1 / 2 j j 1 j 2 1 σ ≅ − ν ( ) c c (1 ) modified flux 2 L-W scheme generates numerical oscillation. Another scheme added non-linear flux is introduced. 1 ~ = + 1 ν − − n f c u [ B ( )( u u )] + + + j 1 2 / j j 1 2 / j 1 j 2 flux limiter function When B select well, monotone scheme, which is almost higher order and becomes 1st order near discontinuity, can be created ・・Modified flux method proposed by Harten 17 17
What is monotone scheme? What is monotone scheme? u monotone solution non-monotone solution x Monotone imply a property to become 1st order near discontinuity 18 18
TVD(Total Variation Diminishing) Method Variation Diminishing) Method TVD(Total Total variation in space at a time defines as follows: ∑ n ) = − u j n n TV ( u u j + 1 : Total variation j TV stability:Total variation decrease with time (diminishing)・・TVD condition n + 1 ) ≤ TV ( u n ) TV ( u Scheme which is satisfied this condition is called TVD scheme. 19 19
TVD(Total Variation Diminishing) Method Variation Diminishing) Method TVD(Total TVD scheme means to preserve monotone profile n + Δ t n + 1 = u j − + ( u j + 1 − u j ) + c j − 1/ 2 ( u j − u j − 1 )] u j Δ x [ c j + 1 / 2 In order to preserve monotone, CFL condition Δ t − + − + > 0, > 0, − c j − 1/ 2 ) ≤ 1 Δ x ( c j + 1/ 2 c j + 1/ 2 c j − 1/ 2 For example, 1st order upwind difference is TVD scheme. monotone TVD scheme 20 20
Recommend
More recommend