Simulation of Turbulent Flows With Strong Shocks

1
Simulation of Turbulent Flows With Strong Shocks

• Bruce Fryxell
• Los Alamos National Laboratory

Collaborators Suresh Menon, Franklin Genin,
Jason Hackl (GIT)
Turbulent Mixing and Beyond Trieste, Italy
August, 2007
2
Objectives and Challenges

• Objective
• To develop a multi-phase (gas, liquid, solid) LES
  solver to study detonation and neutralization of
  chemical/biological agents in ground-fixed
  structures
• Challenges
• Simulation model must capture physics of UNSTEADY
  shock-turbulence-chemistry interactions without
  excessive dissipation
• Shock-induced vaporization and/or ignition of
  dense particles must be included
• Multi-component droplets and multi-species
  finite-rate kinetics
• Complex flow-structure interactions in 3D
• Approach
• Expand capabilities of an existing LES solver
  LESLIE3D to address unsteady shock-turbulence
  physics

3
Baseline Code LESLIE3D

• Well-established DNS/LES code developed at
  Georgia Tech by Suresh Menon and collaborators
  primarily for aerospace and combustion problems
• Extensively verified and validated
• Optimized to run efficiently on a variety of
  massively-parallel computers
• Gas dynamics
• Smooth-flow finite-volume (SF-FV) Navier-Stokes
  solver accurate to second order in time and
  fourth order in space

4
Euler Equations
In one-dimension
Xl is the specific density of the lth fluid
5
Smooth-flow Solver (SF-FV)
Predictor-corrector second-order time
differencing
For second-order in space
For fourth-order in space
Alternating between upwind and downwind fluxes
results in a central difference method
6
LESLIE3D

• Large Eddy Simulation using the ksgs approach
• One additional equation solved for the subgrid
  kinetic energy
• Localized dynamic model (LDKM) for subgrid
  closure of the unresolved momentum and energy
  fluxes that contains no adjustable parameters
• Linear Eddy Mixing (LEM)
• Grid-within-grid approach for reaction-diffusion
  and scalar mixing at the subgrid scale
• Grid-free Lagrangian tracking of liquid and solid
  particles including vaporization and ignition of
  particles
• Finite-rate reaction kinetics
• Thermally perfect and real gas equations of state
• Capable of simulations in complex geometries

7
Algorithm for Strong Shock Propagation

• Baseline gas dynamics solver in LESLIE3D is
  accurate for smooth flows, but cannot treat
  shocks and contact discontinuities without
  generating unphysical oscillations
• Incorporated accurate shock capturing capability
  using the Piecewise-Parabolic Method (PPM) of
  Colella and Woodward
• PPM solver is coupled to all of the other physics
  modules in the baseline code

8
Hybrid Method

• Simulation of flows containing both shocks and
  turbulence presents significant challenges for
  numerical methods
• SF-FV solver can not be used near discontinuities
  without generating unphysical oscillations
• Shock-capturing method could be used to calculate
  entire flow, but the upwind dissipation needed to
  stabilize shocks produces an unphysical rate of
  decay of turbulent features in the flow
• Shock capturing method is much more expensive
  than SF-FV
• Coupling of subgrid models to SF-FV has been
  extensively validated and its properties are well
  understood

9
Hybrid Method

• Solution use a Hybrid Method
• Calculate regions of the flow near
  discontinuities with the shock-capturing method
• Compute remainder of the flow with a high-order
  central difference method (SF-FV)
• Retains the best features of both methods
• Can calculate more accurate answers than could be
  obtained with either method by itself
• More efficient than using the shock-capturing
  method for the entire flow

10
Hybrid Method

• Define a smoothness parameter according to

• Q can be any variable. Here we use both density
  and pressure and take the maximum of the two.
• To prevent triggering on numerical noise, we set
  Si to 0 if either the numerator or denominator
  divided by Qi is less than 0.05.
• For multidimensional flows, the maximum value of
  Si for each coordinate direction is used. Cross
  derivatives are not considered.

11
Hybrid Method

• In any cell in which Si gt 0.35, PPM is used to
  construct the fluxes through each face of that
  cell
• In any cell bordering a cell in which the PPM
  fluxes are used, the average of the PPM and SF-FV
  fluxes are used
• In the rest of the grid, the SF-FV fluxes are
  used
• It should be possible to combine any two methods
  using this technique, although tuning of the
  dimensionless parameters might be required

12
Hybrid Method

• Shu-Osher test problem
• Planar shock propagating through a sinusoidal
  density field
• Tests codes ability to accurately compute shocks
  and short-wavelength smooth variations in the
  same calculation
• Ideal test problem for hybrid code
• Particularly relevant to the problem of
  shock-induced turbulence

13
Initial Conditions
14
Shu-Osher Test Problem
Solid line reference solution Blue SF-FV
fluxes Red PPM fluxes Green Average fluxes
Results for the Shu-Osher problem obtained using
the hybrid method with 400 points at t 1.872
15
Shu-Osher Test Problem
Solid line reference solution Blue SF-FV
fluxes Red PPM fluxes Green Average fluxes
Results for the Shu-Osher problem obtained using
the hybrid method with 800 points at t 1.872
16
Shu-Osher Test Problem
Solid line reference solution Blue SF-FV
fluxes Red PPM fluxes Green Average fluxes
Results for the Shu-Osher problem obtained using
the hybrid method with 1600 points at t 1.872
17
Shu-Osher Test Problem
Solid line PPM 800 pts Blue
Hybrid 800 pts
Closeup of the region immediately behind the
shock in the Shu-Osher problem on an 800 point
grid.
18
Number of Grid Points Where Each Flux is Used
Shu-Osher Test Problem
19
Richtmyer-Meshkov Instability

• Excellent problem for code validation
• Involves both shocks and turbulence for
  sufficiently high Reynolds numbers
• Linear and non-linear analytic theories
• Experimental results
• Code-code comparisons
• Both two and three dimensions
• Single mode and multi mode

20
Initial Conditions
0.04 cm
He
H2
0.01 cm
Shock
Interface
rH 0.001 g cm-3
rHe 0.015 g cm-3
Pshock 300 P0

• Multimode perturbation
• Primary mode has an amplitude of 0.001 cm
• Secondary mode is antisymmetric with 1/5 the
  wavelength and 1/10 the amplitude of the
  primary mode

21
Two Dimensions Multimode - Inviscid
1024 x 256 grid
22
Two Dimensions Multimode -Viscous
1024 x 256 grid
23
Percent of Grid Points Where Each Flux is Used
Two Dimensions Multimode - Inviscid
PPM Fluxes
Averaged Fluxes
SM-FV Fluxes
256 x 64
58
18
24
512 x 128
43
15
43
1024 x 256
23
12
65
24
Two Dimensions Multimode -Viscous
Percent of Grid Points Where Each Flux is Used
PPM Fluxes
Averaged Fluxes
SM-FV Fluxes
256 x 64
46
15
39
512 x 128
32
12
56
1024 x 256
11
10
79
25
Two Dimensions Multimode - Inviscid
Density Gradient
Red PPM Blue SF-FV
26
Three Dimensions Single Mode
Temperature
2D
480 x 120 x 120 grid
Temperature
Red PPM fluxes Blue SF-FV fluxes Green
Averaged Fluxes
3D
In three dimensions the structure of the
instability is significantly different. The
finger is broader and the growth rate is
substantially larger.
27
Three Dimensions Single Mode
Temperature Isosurface
480 x 120 x 120 grid
28
Three Dimensions Single Mode
Temperature Isosurface
480 x 120 x 120 grid
29
Three Dimensions Single Mode
Density Isosurface
480 x 120 x 120
30
Shock - Turbulence Interaction

• Flow of isotropic turbulence through a planar
  shock.
• Level of turbulence increases across shock.
• Want to compare results to DNS calculation of
  Mahesh, Lele, and Moin (1997).
• DNS calculation used 6th order ENO for the shock
  and a 6th order Padé scheme for the turbulent
  flow.
• DNS calculation used 231 x 81 x 81 grid.
• LES calculation uses 62 x 32 x 32 grid.

31
Problem Initialization

• Begin with simulation of decaying turbulence in a
  cubic box of dimensions where k0 is
  the most energetic wave number
• Grid resolution is 32 x 32 x 32
• Initial energy spectrum is given by

32
Problem Initialization

• Initial turbulent Mach number Mt 0.22 decayed
  to a final value of 0.14
• A mean velocity is then added to the turbulent
  flow, and this is used as an inflow boundary
  condition
• Characteristic boundary conditions were used
  downstream
• Periodic boundary conditions were used at side
  boundaries
• Shock has a Mach number of 1.9
• Simulation was allowed to settle down for two
  flow times
• Turbulent statistics were then collected for
  another two flow times

33
Shock-Turbulence Interaction

• 62x32x32 grid points
• High clustering in the near shock region
• Subsonic characteristic outflow conditions

34
Shock-Turbulence Interaction
Subgrid kinetic energy
35
Shock-Turbulence InteractionLongitudinal
Reynolds Stress
36
Shock-Turbulence Interaction Longitudinal
Reynolds Stress
37
Preliminary Agent Defeat Simulations
38
Explosion in Complex Geometry

• Room size 20 m x 15 m x 5 m
• Three objects inserted into room
• Cube 2.5 m on each side located near the back
  and left walls and oriented at a 45o angle
• Wall A 5 m wide, 4 m high, 0.25 m thick located
  near the back and right walls and oriented at a
  30o angle
• Wall B 5 m wide, 4 m high, 0.25 m thick located
  near the front wall and oriented parallel to the
  front wall
• A door 1 m wide and 2 m high is placed in the
  center of the front wall

39
Explosion in Complex Geometry

• Explosion initiated using the analytic solution
  for a Sedov-Taylor blast wave
• Explosion initiated using 2 kg of explosive
• Explosion energy of 4.66 x 106 J/Kg
• Initial diameter of explosion 1.0 m
• Simulation performed on a 320 x 240 x 80 uniform
  Cartesian grid
• Internal objects are treated using the method of
  embedded boundaries of Pember et al.

40
Explosion in Complex Geometry
Pressure 2 m above floor t 3.9 ms
41
Explosion in Complex Geometry
Pressure 2 m above floor t 8.8 ms
42
Explosion in Complex Geometry
Pressure 2 m above floor t 14.9 ms
43
Explosion in Complex Geometry
Pressure 2 m above floor t 21.6 ms