Shock wave - boundary layer interaction in a long shock tube

The numerical simulation of influence of viscous effects on flow inside shock tubes are described. We considered three shock tubes with different length-to-height ratio: 1 to 1 (Daru & Tenaud [6] case), 8 to 1 (extended Daru & Tenaud case) and 116 to 1 (HAST-1 [10-11] tube). Some details of different stages of the shock wave-boundary layer interaction are discussed.


Introduction
Shock tubes are a convenient tool for experimental study of fast-flowing high-speed flows. Interest in them arose long ago and proceeds to present day [1,2]. This is due, first of all, to the complexity and high cost of realization of full-scale tests in high-speed flows. Recently, special attention is spared to carrying out high-precision experiments on the shock tubes for verification and validation of CFD computer codes intended for the numerical simulation of hypersonic flows [3].To improve the accuracy of the experimental results obtained, a clear understanding of the processes occurring in the shock tubes is needed. First of all, this is the influence of viscous effects on the propagation of gasdynamic discontinuities inside the shock tubes. This influence was studied experimentally [3,4] and theoretically [5]. Theoretical studies were conducted on the basis of approximate models such as quasi one-dimensional gas-dynamics and boundary layer equations.
At the present time, the development of computer technologies makes it possible to numerically study the processes inside the shock tubes on the basis of the complete Navier-Stokes equations on very detailed grids and to obtain a detailed structure of the flow in the shock tubes. Most of these numerical studies of the flow in shock tubes concentrate on the investigation of the wave structure that occurs when a shock wave is reflected from a wall. In [6] one such variant was proposed as a benchmark solution and was studied in detail on fine grids for several Reynolds numbers (200, 1000, 5000, etc.) [69]. This variant is characterized by a very small length of the shock tube, which makes it impossible to separate the interaction of various gas-dynamical discontinuities with the boundary layer and with each other.
Real shock tubes are characterized by an appreciably larger length in comparison with a diameter. For example, in the case of HAST-1 [1011] the total length of the high and low pressures chambers is 9.5 m, but the diameter is 0.08 m. This leads to the facts that: 1). The time of propagation of the shock wave from the diaphragm to the shock tube end is large enough for a noticeable influence of viscous effects on the flow inside the shock wave.

Numerical code verification
For our numerical code verification we present some results of numerical simulation of the 2-D viscous shock tube problem. This problem was studied by Daru and Tenaud [6] and Sjögreen and Yee [7]. The 3D version of the problem was proposed in [8]. In the problem a tiny, very short shock tube with square cross section is considered. In 2D the shock tube is a 2D box 0≤x, y≤1. A membrane located at x = 0.5 separates two different states of the gas. At time zero the membrane is removed and wave interaction occurs. The initial state is given as follows: (1) γ = 1.4 and the Prandtl number is 0.73. The two Reynolds numbers considered are 200 and 1000. The viscosity is assumed to be constant and independent of temperature, so Sutherland's law is not used. Figure 1 shows the numerical solution of the problem obtained on a 2000×1000 mesh for Re = 200. Figure 1 contains numerical Shlirien image (the density gradient countours). A comparison with the benchmark solution [69] shows that in this case the key features of the flow dynamics generated by the interaction of the shock wave and the boundary layer are captured well.
The next, we consider the case Re = 1000. For such Reynolds number, vortical structures become more complex. Small vortices are generated just downstream of the boundary layer separation, and shocklets are present at many locations inside the large vortices. Figure 2 presents density gradient magnitude contours at non-dimensional time t end = 1 predicted on mesh 3000×1500. Presented results   Obtained solutions of these two problems show that the flow in shock tube has complex 2-D structure due to shock wave interaction with boundary layers. Usual 1-D numerical simulation of the flow can't reproduce all flow details and therefore has limited applicability.
Comparing these figures, it can be noted that an increase in the Reynolds number leads to a noticeable complication in the flow structure. A large number of small vortex structures appear, for the resolution of which it is necessary to refine of the computational grid in both longitudinal and transverse directions. It is an important reason why the Daru and Tenaud benchmark problem is usually solved on grids with square cells. Therefore the need for a strong refinement limits the possibility of numerical simulation of the flow in shock tubes of large length with high accuracy.

Shock wave propagation in a tube
Let us consider the flow inside the shock tube in more detail. Usually, the rupture of the diaphragm occurs non-uniformly and its further "opening" cannot be complete. This introduces additional perturbations in the flows within shock tube at the initial stage. A very detailed description of the process of opening a diaphragm in a shock tube can be found in [18]. However, at a large distance from the diaphragm, all perturbations decay fairly quickly and the front of the shock wave becomes plane. Therefore, the structure of the flow behind the shock wave is determined by the perturbations arising on the walls of the shock tube.
With sensitive tuning of the interference-shadow instruments in experiment, one can observe a characteristic "net" structure behind the shock wave, which arises as a result of the addition of many random perturbations from the roughness on the wall [17]. Contact surfaces that are part of this structure form vortices which can generate turbulence in the flow behind the shock wave. In our calculations, given below, we shall assume that the Reynolds number is sufficiently small and there is no transition to turbulence. The second important assumption, which confirms the laminar nature of the flow, is the smoothness of the shock tube walls, so there are no additional "micro-perturbations" that can initiate the transition.
A similar assumption is made about the initial rupture of the diaphragm: the diaphragm rupture is modeled by an instant complete "disappearance" of the diaphragm, which does not introduce additional perturbations into the flow.
If viscous effects are not taken into account, then in shock tube, after a diaphragm rupture, a flow close to the one-dimensional solution of the Riemann problem is realized [19]. As a rule, the ratio of pressures in the chambers of high pressure (to the left of the diaphragm) and low pressure (to the right of the diaphragm) leads to the formation of the following configuration of discontinuities: the shock wave moves to the right, followed by a contact discontinuity separating the gases in the chambers of high and low pressures. The main difference in pressure is compensated for by a large fan of rarefaction waves, which propagates in both directions, but mainly to the left. The initial state in the Daru and Tenaud problem [6] gives the same flow configuration after a diaphragm rupture.
Part of the rarefaction fan moving to the left reaches the wall, is reflected from it and begins to move in the direction of the contact discontinuity and shock wave. Velocity of the reflected part of the rarefaction fan is higher than velocity of the right boundary of the original fan. As result, under a certain initial conditions, the reflected fan can "catch up" with the right boundary of the original fan and begin to influence the flow structure near the "working" right end of the shock tube.
It will be convenient to consider the following stages of the flow in the shock tube separately: 1). Motion of the initial configuration of discontinuities along the shock tube before the interaction of the incidental shock with the wall.
2). Bifurcation of reflected shock wave before the interaction of the reflected shock with the contact discontinuity.
3). The interaction of the reflected-shock bifurcation with the contact discontinuity until the arrival of the right rarefaction wave from the fan.
For numerical simulation a shock tube with a length-to-height ratio of 8 to 1 was considered. The high-pressure chamber occupied the left part of the shock tube [0,0.5) x  , and the low-pressure chamber occupied the right one (0. 5,8] x  . To construct a computational grid, we used the grid step . So one-dimensional calculations were performed on a grid of 8000 cells, and two-dimensional calculations were performed on a grid of 8000 × 500. The initial state of the gas in the shock tube corresponds to the Daru and Tenaud problem (1).
If the flow in a shock tube is considered in the 1D formulation or the 2D inviscid formulation, then the flow parameters behind the shock front must be constant up to the contact surface, the speed of which is equal to the gas velocity behind the shock wave. This is clearly seen in figure 3, where the left part of the figure shows the density along the shock tube axis at some time, and the right part shows the pressure. It is seen that in one-dimensional calculations the density and pressure between the shock waves and the contact discontinuity are constant. In reality, the flow parameters behind the shock wave do not remain constant, as is clearly seen in the pressure plot in figure 3 in the case of a two-dimensional viscous calculation. This is due to the fact that the boundary layer on the walls of the shock tube grows over the entire interval from the shock wave to the contact discontinuity (see figure 4). From figure 4 it can be easily seen that the growth of the boundary layer leads to a change in the "effective cross-section" of the shock tube and, consequently, to a change in the parameters in the inviscid core of the flow [20]. The flow velocity increases smoothly from the shock wave to the contact surface [21], the temperature and pressure increase. These changes can be estimated from the approximate theory [5].
From figure 4 it can be also seen that in the viscous 2D case the shock speed decreases, and the contact surface speed increases in comparison with the 1D case when these speeds are constant.
We shall note that in long shock tubes, the boundary layers on the side walls can merge. In this case, one-dimensional calculations are not very suitable for estimating parameters in shock tubes.
It should also be noted that taking into account viscous effects in the longitudinal direction in the 1D approximation in this case has very little effect on the flow in the shock tube, as is clearly seen in figure 3.

Bifurcation of reflected shock wave
When the shock wave is reflected from a rigid wall in a shock tube, it propagates along an inhomogeneous flow. Mainly this flow consists of two layers: the inviscid core and the boundary layer formed after the incident shock wave. Griffiths predicted [22], and Mark confirmed experimentally [23] that flow regimes with the bifurcation of reflected shock can arise in the viscous flow, and moreover the characteristic distance from the wall to the splitting point of the reflected shock is greater than the thickness of the boundary layer. The bifurcation of reflected shock in a shock tube can be obtained numerically too. Some complicated variants of the bifurcation of reflected shock in the Daru and Tenaud shock tube are shown in figures 1 and 2. Schematically, a lambda-shape like bifurcated shock wave pattern is shown in figure 5. After reflection the stagnation pressure of the boundary-layer becomes lower than the pressure behind the reflected shock wave. As result the fluid cannot pass under the reflected shock and a separated flow region appears. This separated region propagates upstream with the reflected shock ( figure 5).
The problem of calculating the wave structure arising during bifurcation is not completely solved. The main problem is the obtaining of a convergent solution when the grid is refined [9] in the case of large Reynolds numbers. This was one of the reasons for choosing Re = 200 in this numerical simulation. The temperature field and the Mach number contour lines are shown in figure 6 after the shock reflection. The lambda-shape like bifurcated shock structure is clearly visible, which completely corresponds to figure 5. In addition, in figure 6 on the left, the contact surface is also visible, moving towards the reflected shock wave.  It should be noted that there is a slight decreasing of pressure near the end of the shock tube in the two-dimensional case in comparison with the one-dimensional calculation.    The speed of the triple point movement is quite high and very soon the triple point reaches the tube axis (see figure 9). It can be seen the bifurcation shock junction. The flow near the shock tube end becomes two-dimensional.

Numerical simulation of the flow in a real shock tube
In real shock tubes such as described in [1011] the length is considerably large than the tube height. The total length of driver and driven parts of HAST1 [1011] equals to 9.32 m and the diameter of this shock tube is 0.08 m. It means that the length-to-diameter ratio is 116.5. With such large ratio the possibility exists that the boundary layers on the tube walls are closed. This may affect the structure of the shock wave and the dynamics of its propagation.
Such large length-to-height ratio introduces additional computational difficulties. It doesn't allow to generate a computational grid with the same resolution as for Daru and Tenaud problem above. We used a 5000×200 uniform mesh in our numerical simulation. The comparison this solution with the solution on a 4000×150 mesh shows that all base features of the flow are very similar on the both meshes. In the driver tube we set the pressure 20 atm and the pressure in the driven tube was 10 -3 atm. The temperature in both tubes was 300˚K. Figure 10 shows the density gradient magnitude contours at different times. It should be noted that in this figure essentially different length scales along different axes are used. On the top part of the figure the density gradient contours at a time shortly before the shock wave interaction with the wall are shown. It can be seen considerable growth of the boundary layer which takes up a third of the tube radius. The remaining three pictures correspond to three different times after the shock wave reflection from the wall. aren't resolved in the presented solution. Insufficiency of the grid resolution becomes even more pronounced when reflected from the left wall rarefaction fan reaches the right end of the tube and begins to interact with the reflected shock wave.

Conclusion
The numerical simulation of influence of viscous effects on flow inside shock tubes are described. We considered three shock tubes with different length-to-height ratio: Daru & Tenaud [6] case (1 to 1 ratio), extended Daru & Tenaud case (8 to 1 ratio) and HAST-1 tube [1011] (116 to 1 ratio). The Daru & Tenaud case was used as a benchmark problem. Good agreement between the obtained results and well-known solutions was demonstrated. The extended Daru & Tenaud case allowed us to separate different stages of the unsteady flow in shock tube and compare results of 1D and 2D calculations.
The numerical simulation of the flow in HAST-1 shows the significant influence of viscous effects on the flow.