Open Journal of Fluid Dynamics
Vol.08 No.04(2018), Article ID:88433,12 pages

Numerical Simulations of Shock-Shock Interactions

Roberto Paoli1,2

1Department of Mechanical and Industrial Engineering, University of Illinois at Chicago, Chicago, IL, USA

2Argonne National Laboratory, Computational Science Division, Argonne, IL, USA

Copyright © 2018 by author and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

Received: August 21, 2018; Accepted: November 10, 2018; Published: November 13, 2018


We present a computational study of Edney type shock-shock interactions between an oblique shock and a bow shock in front of a circular cylinder. Both cold and hot hypersonic test cases were considered and the results revealed the main features characterizing the type III and type IV of interaction, without and with the jet impinging on the cylinder. The simulations also predicted qualitatively the complex flow structure observed in the experiments and the pressure and heat peak on the surface of the cylinder.


Shock Waves, Edney Type Interactions, Non-Equilibrium Flows

1. Introduction

In the present work we analyze the effects of thermochemical relaxation phenomena on the interference between a cylindrical bow shock and a plane oblique shock. This problem has been extensively studied since the 60’s [1] [2] because of its importance in many aerospace applications involving complex shock-shock interactions such as fin-fuselage interaction in re-entry vehicles, launcher-booster interaction as well as cowl-compression ramp shock interaction in engine inlets.

According to Edney [1] , for a given incident shock angle, the shock wave configuration consists of six possible interference patterns, depending on the position of the impinging shock with respect to the bow shock. In particular, it is known that Type III and IV configurations are the most critical, being characterized by a severe increase in the dynamic and thermal loads on the body surface (see the sketch in Figure 1).

In both cases the shock pattern is characterized by the occurrence of two triple

Figure 1. Sketch of the shock-shock interaction according to Edney [1] . Adapted from Sanderson [9] .

points caused by the interaction of the incident shock with the bow shock. It should be remarked that in a Type-III interference the amplification of the surface pressure and heat transfer is caused by the attached shear layer while in a Type IV the amplification is mainly due to the impingement of the supersonic jet originating at the lower triple point and terminating with a “nearly normal” jet shock in front of the body.

For cold hypersonic flows many works (both numerical and experimental) can be found in the literature. Wieting and Holden [3] have experimentally studied shock-shock interaction at M = 6 and 8 at various Reynolds numbers. Those authors have observed an extremum in the Type-IV interaction both in surface pressure and heat flux peaks. They have also pointed out that, for high Reynolds numbers, the heat flux exhibits an additional increase due to the jet turbulence. Borovoy et al. [4] have reported an experimental study of shock/shock interactions at M = 6 and have observed the occurrence of a double peak in the heat flux as a function of the shock impingement location due to a change of the jet structure near the transition from a Type-III to Type-IV interaction. In addition, they have also observed a weak dependence of the thermal load upon the Reynolds number. Grasso et al. [5] have conducted a theoretical and experimental investigation of the heat transfer rates both in absence and in the presence of shock interaction. In the former case, the experiments showed that acoustic phenomena may amplify laminar-to-turbulence transition and consequently increase the peak heating. In the presence of shock/shock interaction, they observed a displacement of the peak loads (pressure and heat) locations; in addition, from comparison between experimental results and theoretical predictions, they concluded that the correlations for pressure are most sensitive to the deformation of the bow shock subsequent to shock impingement.

Considering now real gas effects, it is known that, when the free-stream specific kinetic energy is of the same order of magnitude of the dissociation energy, non-equilibrium phenomena must be taken into account [6] [7] [8] [9] [10] . The main effect of thermochemical relaxation is the increase in the density ratio across the shock waves thus causing a reduction of the length scales (in particular, the bow shock stand-off distance and, for a Type-IV interaction, the jet width and jet shock stand-off distance). The effects of relaxation behind strong shock waves in a Type-IV interaction have been studied both theoretically and experimentally by Sanderson [9] who observed that the difference between the equilibrium and the frozen heat transfer amplifications decreases as the Mach number increases. Another critical aspect is that, for some particular flow conditions (i.e. shock strength and Mach number), a Type-IV interaction exhibits some unsteadiness due to the vortex shedding that takes place in the terminating part of the supersonic jet, and whose frequency depends on the characteristic scales of the flow-field such as stand-off distance, length and width of the jet [10] [11] [12] . Indeed, relaxation phenomena contribute to increase the vortex shedding frequency, thus strongly affecting the unsteady character of the flow. In addition, depending upon the value of the free-stream stagnation enthalpy, the recombination rate may be sufficiently high for the boundary layer to be in equilibrium thus contributing to increase the peak heat transfer.

The motivation of the present work is to assess the influence of non-equilibrium on Type-III and Type-IV shock/shock interactions using numerical modeling In particular, we are interested on the effects of non-equilibrium on the flow scales and unsteadiness, and on the dynamic and thermal loads. This represents a first step where we focus on the complex shock pattern in the case of ideal gas simulations and present preliminary results of the simulations for thermochemical nonequlibrium. The approach relies on a finite volume method for the solution of the conservation laws for a mixture of gases in non-equilibrium. In addition, in order to understand the effects of vibrational relaxation and chemical reactions we have also performed ideal gas simulations.

In this section we formulate the two dimensional conservation equations for a multispecies non-ionizing gas in thermochemical non-equilibrium, following the approach developed in Grasso and Capano [13] and Grasso and Paoli [14] . If one introduces a single vibrational temperature to characterize the internal energy mode of all diatomic species, the conservation equations can be cast in vector form as

t S w d S + S f n d l = S w ˙ d S (1)


w = [ ρ k , ρ u , ρ v , ρ E , ρ e v ] T

f = f E f V

f E = [ ρ k u , ρ u u + p I , ρ u H , ρ u e v ] T (2)

f V = [ J k , σ , σ u J q , J q v ] T

w ˙ = [ ρ ˙ s , 0 , 0 , 0 , ρ e ˙ v ] T

and the mass flux of species k ( J k ), the internal ( J q ) and vibration ( J q v ) energy diffusion fluxes (due to conduction and species diffusion) and the stress tensor σ are

J k = ρ D k Y k (3)

J q = j q t r + j q v + k h k J k = ( η t + η r ) T η v T v + k h k J k (4)

J q v = j q v + k h k v J k = η v T v + k h k v J k (5)

σ = μ ( u + u T ) 2 3 μ u I (6)

where k indicates diatomic species and Y k = ρ k / ρ is the species mass fraction, while the transport coefficients ( D k , η’s and μ ) are defined according to Chapmann-Enskog theory [13] .

For atomic species the internal energy ( e k ) accounts only for the translational ( e k t ) and the enthalpy of formation ( Δ h k 0 ) contributions

e k = e k t + Δ h k o = 3 2 R k T + Δ h k o (7)

while for diatomic species the internal energy also accounts for rotation ( e k r ) and vibration ( e k v ) contributions and

e k = e k t + e k r + e k v + Δ h k o = 3 2 R k T + 1 2 R k T + R k θ k v exp ( θ k v / T v ) 1 + Δ h k o (8)

where θ k v is the characteristic vibrational temperature. The vibrational energy ( e v ) and the total enthalpy (H) of the mixture are

e v = k Y k e k , H = E + p ρ = k Y k e k + u u 2 + p ρ (9)

The species production ( ρ ˙ k ) due to finite rate chemistry has been modeled by means of Park’s reaction mechanism whereby we have introduced a rate controlling temperature for the dissociation reactions defined as T d = T T v so as to account for the coupling between vibration and dissociation. The vibrational energy source ( e ˙ v ) is the sum of the the T-V energy exchange contribution ( S T V ) and to the energy removal due to vibration-dissociation coupling ( S V D ) that have been modeled as discussed in Grasso and Capano [13] according to

S T V = k ρ k e k v * e k v τ k v , S V D = k ρ ˙ k e k v (10)

where e k v * = R k θ k v / ( exp ( θ k v / T ) 1 ) is the equilibrium vibrational energy of molecular species, and τ k v is the relaxation time obtained from the Millikan-White relation [13] .

2. Numerics

The basic numerical algorithm relies on a cell-centered finite-volume “patched” subdomain formulation and the governing equations are cast in the following fully discretized form:

S i , j d w i , j d t + β = 1 4 ( f ^ E f ^ V ) β n β Δ l β = S i , j w ˙ i , j (11)

where Δ l β and S i , j are, respectively, the size of face β and the cell area, while f ^ E and f ^ V represent, respectively, the numerical inviscid and viscous fluxes. The former is based on a second order upwind biased total variation diminishing scheme that accounts for non-equilibrium phenomena, following the general methodology described in Grasso and Capano [13] .

For example, the numerical inviscid flux discretization at the generic cell face ( i + 1 2 , j ) , is cast in the following form:

( f ^ E n ) i + 1 2 , j = 1 2 [ ( f E n ) i , j + ( f E n ) i + 1 , j + R i + 1 2 , j Φ i + 1 2 , j ] (12)

where Φ is obtained by characteristic decomposition in the direction normal to the cell face and the use of a minmod slope limiter. The right eigenvector matrix is defined as

R = [ δ k q 0 Y k Y k Y k u c n y 0 u + c n x u c n x v c n x 0 v + c n x v c n y u u / 2 χ q / K c ( u b ) 1 H + c n x H c n y 0 0 1 e v e v ] (13)

where b n = 0 and the pressure derivatives K and χ q are defined as

K = k Y k R k k Y k c v k t ; χ q = R q T K e q t

and c v k t and c 2 are, respectively, the species translational-rotational specific heat and the frozen speed of sound that are defined as

c v k t = e k t T ; c 2 = q χ q Y q + K ( h e v )

The values at cell interfaces are calculated by using a generalization of Roe’s averaging [15] to allow for thermal and chemical non-equilibrium (see Grasso and Capano [13] for details).

In addition, in order to avoid the carbuncle phenomenon which takes place in all Roe based schemes [16] [17] we adopted the multidimensional dissipation correction to the original dissipation function, following the approach proposed by Morano et al. [16] .

Finally, the time integration is performed by a three-stage Runge-Kutta algorithm where the source terms are treated by a point implicit algorithm by introducing a precondition matrix which is related to the partial Jacobian of the source term [13] .

3. Results

The objective of the present paper is to analyze Type-III and Type-IV shock-shock interactions with the main intent of studying the evolution of the shock wave pattern consequent to the variation of the shock impingement location (steady simulations). We assumed the free-stream conditions of the experimental work of Grasso et al., [5] which, indeed, produced a steady flow configuration. Due to the low pressure and temperature ( p = 5.9 Pa , T = 52 K ), no relaxation effects take place, and we have only performed ideal gas simulations.

For the problem under investigation the predictions of the peak (thermal and dynamic) loads and their locations are mainly affected by the mesh spacing along the body (ξ), provided the minimum grid spacing in the normal-to-wall direction (η) is of the order of microns. We have then performed a grid independence study by changing the mesh spacing in the ξ-direction. In particular we have considered three grids consisting of 144 × 160 , 192 × 160 and 294 × 160 total number of cells (having a minimum normal-to-wall spacing of 3.34 × 10−6 m). The results of the analysis (not reported) showed that grid independency is achieved with the 192 × 160 grid, which has then been selected for all simulations.

In order to span the Type-III and Type-IV shock wave interference patterns, we have followed the same quasi-steady procedure adopted in the experiment by moving the incident shock generator in the vertical direction so as to change the location of the impingement location on the cylinder surface (for a fixed shock angle β = 23.57 deg ). In particular, the shock generator is displaced vertically once a steady (or pseudo steady) flow configuration (for the given shock impingement location) is attained.

In Figure 2 we report the numerical Schlieren plots of three shock interference patterns corresponding to three different shock impingement location on the cylinder (identified in terms of the angular position on the cylinder, ϕ i ). In particular, the figures shows a Type-III configuration ( ϕ i = 9.75 deg ) where the attached shock, due to the reflection of the shear layer on the surface, is clearly visible; a Type-IV configuration ( ϕ i = 14.08 deg ), characterized by the detachment of the reflected shock and the consequent occurrence of a normal shock in front of the body. Finally, the figure shows a Type-IVa configuration which is characterized by a supersonic jet grazing the surface rather than impinging on the cylinder. The effects of the transition from Type-III to IV shock pattern on density flow-field are also shown in more detail in Figure 3 where we report the numerical interferogram around the impingement region and compare with the experimental interferogram for a similar configuration. In particular, one observes the rapid increase in the maximum (normalized) density level as the flow passes from Type-III and IV shock structure, as an increase of the impingement of the supersonic jet on the body surface. The evolution of the Type-IV structure in the neighborhood of the transition limit is shown in Figure 4 where we report the iso-density flow field around the terminating jet shock. In order to clearly identify the boundary of the supersonic jet region as well as the relative position between the body and the jet shock, the maximum density level displayed has been set equal to the maximum value attained inside the jet, and

Figure 2. Numerical Schlieren of various interference patterns ( M = 9.95 , β = 23.57 deg ). Top left: Type III, ϕ i = 9.75 deg ; top right: Type IV, ϕ i = 14.38 deg ; bottom: Type IVa, no impingement. Axis coordinates normalized by the reference length L, corresponding to the radius of the cylinder.

Figure 3. Numerical interferograms of Type-III and Type-IV configurations ( M = 9.95 , β = 23.57 deg ). Top: Type III, ϕ i = 9.75 deg ; bottom: Type IV, ϕ i = 14.38 deg . The experimental interferograms from Sanderson [9] are also reported for qualitative comparison in the right panels.

Figure 4. Evolution of the shock interference pattern ( M = 9.95 , β = 23.57 deg ) with enlarged view around the impingement location. Normalized density contour plots and, superimposed, M = 1 contour line. Top left: Type III; top right: Type IV and ϕ i = 8.05 deg ; bottom left: Type IV and ϕ i = 0.1 deg ; bottom right: Type IV and ϕ i = 14.38 deg .

in addition, we also report the sonic iso-Mach line. We observe that, as the shock impingement location moves upward, the jet shock inclination with respect to the body decreases. Consequently, the pressure and thermal loads increase as shown in Figure 5, where surface pressure and heat transfer are reported as a function of the angular position on the cylinder. In addition, the interaction of the shock originating at the lower triple point with the upper shear layer causes the occurrence of a lambda shock structure consisting in a reflected shock and a Mach stem as also found by Pandolfi and D’Ambrosio [17] . We argue that this is due to the finite thickness of the shear layer, which is responsible for the curvature of the shock and, ultimately, for the formation of the lamba shock

Figure 5. Surface pressure (left) and heat (right) amplifications (with respect to the values in the absence of interaction).

structure; note that such phenomenon is enhanced as the Mach number increases since the interaction takes place closer to the body where the shear layer is more spread thus leading to higher pressure gradients.

We finally present preliminary results of simulations carried out with air in thermochemical non-equilibrium. The stagnation quantities were taken from the experimental investigation of Carl et al. [10] which provided the following free-stream conditions: u = 5700 m / s , ρ = 2.19 g / m 3 , T = 797 K ; consequently, the free-stream frozen Mach number was M = 10.04 . Ideal gas simulations were finally performed by assuming the same free-stream frozen Mach number. The results for non-equilibrium flow computations are reported in Figure 6 in terms of frozen Mach number and in terms of vibrational temperature and oxygen mass fraction, respectively. The relaxation processes are clearly visible across the strong shock waves, while the weak incident shock remains essentially frozen.

Figure 6. Simulations with real gas effects ( M = 9.95 , β = 23.57 deg ). Top left: frozen Mach number; top right: normalized vibrational temperature; bottom: oxygen mass fraction.

4. Conclusion

In this study we presented and discussed the results of numerical simulations of shock-shock interactions between an oblique shock and a bow shock in front of a circular cylinder. These interactions are categorized as types III and IV according to Edney experimental work. The simulations revealed the main features characterizing the two types of interaction and also predicted the pressure and heat peak on the cylinder. The study was conducted mainly for the cold hypersonic test case of Ref. [5] and preliminary analysis of hot hypersonic test case of Ref. [10] was also presented. This revealed the non-equilibrium phenomena occurring when real-gas effects were activated such as thermodynamic vibrational excitation and oxygen dissociation. In the future, we will analyze the interaction in a fully 3D configuration in the presence of thermochemical non-equilibrium.


The author wishes to thank Francesco Grasso for initiating this work.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

Cite this paper

Paoli, R. (2018) Numerical Simulations of Shock-Shock Interactions. Open Journal of Fluid Dynamics, 8, 392-403.


  1. 1. Edney, B.E. (1968) Anomalous Heat Transfer and Pressure Distributions on Blunt Bodies at Hypersonic Speeds in the Presence of an Impinging Shock. FFA Report, Stockholm.

  2. 2. Hains, F.D. and Keyes, J.W. (1972) Shock Interference Heating in Hypersonic Flows. AIAA Journal, 10, 1441-1446.

  3. 3. Wieting, A.R. and Holden, M.S. (1989) Experimental Shock-Wave Interference Heating on a Cylinder at Mach 6 and 8. AIAA Journal, 27, 1557-1565.

  4. 4. Borovoy, V.Y., Chinilov, A.Y., Gusev, V.N., Struminskaya, I.V., Delery, J. and Chanetz, B. (1997) Interference between a Cylindrical Bow Shock and a Plane Oblique Shock. AIAA Journal, 35, 1721-1728.

  5. 5. Grasso, F., Purpura, C., Delery, J. and Chanetz, B. (2003) Type III and Type IV Shock/Shock Interferences: Theoretical and Experimental Aspects. Aerospace Science and Technology, 7, 93-106.

  6. 6. Fay, J.A. and Riddel, F.R. (1958) Theory of Stagnation Point Heat Transfer in Dissociated Air. Journal of Aeronautical Sciences, 25, 2.

  7. 7. Hornung, H.G. (1972) Non-Equilibrium Dissociating Nitrogen Flows over Spheres and Circular Cylinders. Journal of Fluid Mechanics, 53, 149-176.

  8. 8. Wen, C.Y. and Hornung, H.G. (1995) Non-Equilibrium Dissociating Flow over Spheres. Journal of Fluid Mechanics, 299, 389-405.

  9. 9. Sanderson, S.R. (1995) Shock Wave Interaction in Hypervelocity Flow. Ph.D. Thesis, Caltech.

  10. 10. Carl, M., Hannemann, V. and Eitelberg, G. (1998) Shock/Shock Interaction Experiments in the High Enthalpy Shock Tunnel Goettingen. AIAA Paper 98-0775.

  11. 11. Lind, C.A. and Lewis, M.J. (1996) Computational Analysis of the Unsteady Type IV Shock Interaction of Blunt Body Flows. Journal of Propulsion and Power, 12, 127-133.

  12. 12. Furumoto, G.H., Zhong, X. and Skiba, J.C. (1997) Numerical Studies of Real-Gas Effects on Two-Dimensional Hypersonic Shock-Wave/Boundary-Layer Interaction. Physics of Fluids, 9, 191.

  13. 13. Grasso, F. and Capano, G. (1995) Modeling of Ionizing Hypersonic Flows in Non-Equilibrium. Journal of Spacecraft and Rockets, 32, 217-224.

  14. 14. Grasso, F. and Paoli, R. (2000) Stability of Shock Wave Reflections in Non-Equilibrium Steady Flows and Hysteresis. Physics of Fluids, 12, 3265.

  15. 15. Roe, P.L. (1981) Approximate Riemann Solvers, Parameters Vectors, and Difference Schemes. Journal of Computational Physics, 43, 357-372.

  16. 16. Sanders, R., Morano, E. and Druguet, M.C. (1998) Multidimensional Dissipation for Upwind Schemes: Stability and Applications to Gas Dynamics. Journal of Computational Physics, 145, 511-537.

  17. 17. Pandolfi, L. and D’Ambrosio, D. (2001) Numerical Simulations of Shock Interactions. Proc. of Computational Fluid Dynamics for the 21st Century, Kyoto, Japan, 15-17 July 2001, 322-335.