Simulation of Non-Ideal Effects in Shock-Tubes Gauss Centre for Supercomputing e.V.


Simulation of Non-Ideal Effects in Shock-Tubes

Principal Investigator:
Andreas Kempf

Chair for Fluid Dynamics, University of Duisburg-Essen

Local Project ID:

HPC Platform used:
Hazel Hen of HLRS

Date published:

Scientific work accomplished and results obtained

Non-reactive shocktube simulations Sub Project 1 (SP1)

The aim of the first sub-project was the simulation of non-reactive shocktube experiments for which non-ideal behaviour was observed. This non-ideal behaviour can be caused by shock attenuation [1], or the appearance of bifurcations, due to the interaction of the reflected shock with the boundary layer, which developed behind the incident shock.

When an inert gas, like Argon, fills the driven part of the shocktube, a bifurcation usually does not emerge [2] and the effects of shock attenuation can be investigated exclusively. In most of these setups, an increase of pressure and temperature in time behind the reflected shock occurs, induced by higher temperatures and velocities behind the incident shock upstream.

However, recent experiments showed a decrease of pressure p5 at low pressures.

Test times of the experiments spanned several ms, requiring a very long computational domain. Hence, the simulations were performed in 2 dimensions. The decrease of pressure from the experiments could not be validated yet. Hence, this subproject has not been finished and further simulations will be performed.

Reactive shocktube simulations in 2D Sub Project 2 (SP2)

Simulations of the end-part of a shocktube, housing a stoichiometric H2/O2 mixture have been conducted within this subproject. The numerical setup was set according to the experiments presented in [3]. Depending on the incident Mach-number, mild ignition or strong ignition was observed. The shocktube used for the experiments had a rectangular cross section (1:25 x 1:75 inch). The shorter dimension of the shocktube was chosen for the ”channel-height” of the 2D simulations.

The goal of this sub-project was to investigate the effects of grid resolution, numerical schemes and initial conditions on the ignition delay time Tig in the event of mild ignition and also to prepare larger simulations in 3D. Snapshots of Pseudo Schlieren, superimposed with temperature fields above a threshold value of 1000 K, are presented in fig. 1.

Three different grid resolutions have been tested: 100 μm [a], 50 μm [b] and 25 μm [c]. The corresponding ignition delay times were 170 s [a], 175 s [b] and 177 s [c]. Though huge discrepancies were observed for the three grid resolutions regarding the amount of wave phenomena and shock patterns, the ignition delay times were quite close, suggesting an underlying mechanism responsible for the event of mild ignition, which is not sensitive to grid resolution.

To examine the sensitivity of the mild ignition process on the initial conditions, velocity perturbations were added to to the initial velocity field with a small standard deviation of 1 m/s [e]. Unsurprisingly, the ignition delay time (Tig = 175 μs) [d] did not change, but the exact location of the ignition kernel slightly changed emphasizing the stochastic nature of mild ignition. The sub-project 2 has been finished successfully and the main findings can be used by other groups, who can not afford the computational costs of 3D shocktube simulations or do not have a suitable code for large scale simulations.

Reactive shocktube simulations in 3D, Sub Project 3 (SP3)

The same experimental setup, as in SP2, was used in this subproject. Here, the simulations were carried out in three dimensions to consider realistic turbulence, while being far more expensive with respect to computational resources. The ongoing 3D simulations simulate the complete endwall part of the shocktube.

The main goal of this subproject was to compare the results to those of the previous 2D simulations and to find the main mechanism leading to mild ignition for this specific setup. For this purpose, additional Lagrangian particles were initialised randomly behind the reflected shock, tracking the local state in time.

Figure 2 presents fields of temperature and heat release of a simulation, using a symmetry boundary condition, with a grid resolution of 100 μm. Well visible is the increased heat release in the region of mild ignition on the right of the image. Turbulence near the symmetry boundary condition seems to develop nonphysically. Volume rendered fields of Pseudo Schlieren superimposed with temperature fields, similar to those shown in SP2, are presented in fig. 3. The ignition in the present case happened after 305 μs, hence 3.3 times earlier than in the case of strong ignition, needed for meaningful measurements. Shortly after ignition, the flame front develops into a detonation, as can be seen in the second and third image.

The time- and pressure histories of the Lagrangian particles are illustrated in fig. 4. The time history of the particle, which was located within the ignition kernel at the ignition time, is colored red. For comparison, time histories of particles located near the ignition kernel but further towards the end-wall are presented in black color and particles located in the vicinity of the end-wall are presented in orange color.

It is striking that the initial temperature peak behind the reflected shock increases further upstream. This indicates an acceleration of the reflected shock, since the effect of attenuation at the inlet is not considered in the present simulations.

In order to prove the acceleration of the reflected shock, particle data were used to reconstruct the distance of the reflected shock to the end-wall in time, using the particle location at their maximum temperature. The result is shown in fig. 5. Apparently, the shock accelerates at a certain distance to the end-wall, reaching a velocity offset of roughly 60 m/s, until it decelerates to the original speed. This finding is really interesting, since temperature offsets behind the reflected shock are usually contributed to shock attenuation of the incident shock, as reported for example in [1]. However, the reasons for the acceleration of the reflected shock in the present case need further investigation.

Realization of the project

All the simulations in this project used the in-house code PsiPhi. Throughout the simulations, a cell width of Δ = 100 μm or below was used to have a good representation of boundary layer development. The main runs needed several smaller initial runs. The first run always resembled the typical Riemann problem as initial solution. After the start of the simulation, an expansion wave and a shock developed, while only the shock was of interest for the main run. Hence, at the end of the first simulation, the solution around the shock was stored, in order to be used in the following runs as initial solution. In some cases, a channel flow simulation with periodic inlet/outlet was initialised with the state immediately behind the incident shock, to simulate the correct development of the boundary layer. The mean quantities and the variance of the velocities were calculated for each location in the slice normal to the axial direction and stored as function of time. This allowed to introduce pseudo turbulence at the inlet of the main runs.

Further project specific details are given below.

Non-reactive shocktube simulations SP1

Sub-project SP1 required two main runs in 2D, each utilizing a total of 4096 cores. A 50 by 50
grid was used per Rank, amounting to a total of 10.3 M cells.

Reactive shocktube simulations in 2D SP2

Sub-project 2 required several runs, with a number of cores ranging from 1,440 to 23,040. Accordingly, the total number of cells ranged from 0.6 to 10 M. Problems arose due to the stiffness of the chemistry. A fast, semi-implicit scheme for the time integration of the chemistry was not stable for all the states in the shocktube during the simulation. Sub-time stepping on the other hand required too many computational resources. An implicit solver provided in Sundial libraries finally proofed to be sufficiently fast and stable.

Reactive shocktube simulations in 3D SP3

Several runs have been performed for sub-project SP3, using different grid resolutions, boundary conditions and states behind the reflected shock. A 3D simulation with a grid resolution of Δ = 50 μm utilized 19,800 cores for 48 h (950,400 core-h), if a symmetry boundary condition was applied. The largest simulation, where the computational domain covered the complete end-part of the shocktube at a grid resolution of Δ = 50 mum, ran on 62,424 cores for 24h (1.5 M core-h).


[1] Eric L. Petersen and Ronald K. Hanson, Nonideal effects behind reflected shock waves in a high-pressure shock tube, Shock Waves, 10, no. 6, 405–420, Jan 2001.

[2] Herman Mark, The interaction of a reflected shock wave with the boundary layer in a shock tube, National Advisory Committee for Aeronautics, 1958.

[3] J.W. Meyer and A. K. Oppenheim, On the shock-induced ignition of explosive gases, Symposium (International) on Combustion, 13, no. 1, 1153 – 1164, 1971.

[4] A. Suresh and H. T. Huynh, Accurate Monotonicity-Preserving Schemes with Runge–Kutta Time Stepping, Journal of Computational Physics, 136, no. 1, 83 – 99, 1997.

[5] A. Khokhlov, J. Austin, and A. Knisely, Development of Hot Spots and Ignition Behind Reflected Shocks in 2H2 + O2, Proc. 25th Int. Colloq. Dynamics of Explosions and Reactive Systems. ICDERS, Leeds, UK, Aug. 2015.

Scientific Contact

MSc. Timo Lipkowicz
Chair for fluid dynamics
University of Duisburg-Essen
Carl-Benz-Straße 199, D-47057 Duisburg (Germany)
e-mail: timo.lipkowicz [@]

HLRS Project ID: GCS-snef

July 2019

To top

Tags: University of Duisburg-Essen Computational and Scientific Engineering HLRS