Numerical Investigation of a Shock Accelerated Heavy Gas Cylinder in the Self-Similar Regime

A detailed numerical simulation of a shock accelerated heavy gas (SF6) cylinder surrounded by air gas is presented. It is a simplified configuration of the more general shock-accelerated inhomogeneous flows which occur in a wide variety of astrophysical systems. From the snapshots of the time evolution of the gas cylinder, we find that the evolution of the shock accelerated gas cylinder is in some ways similar to the roll-ups of a vortex sheet for both roll up into a spiral and fall into a self-similar behavior. The systemic and meaningful analyses of the negative circulation, the center of vorticity and the vortex spacing are in a good agreement with results obtained from the prediction of vorticity dynamics. Unlike the mixing zone width in single-mode or multi-mode Richtmyer-Meshkov instability which doesn’t exist, a single power law of time owing to the bubble and spike fronts follow a power law of tθ with different power exponents, the normalized length of the shock accelerated gas cylinder follows a single power law with θ = 0.43 in its self-similar regime obtained from the numerical results.


Introduction
When an impulsive acceleration impinges on the corrugated interface between two fluids of different densities, the instability at the interface will arise due to the deposited vorticity induced by the baroclinic torque produc-tion term 2 p ρ ρ ∇ × ∇ (where p is the pressure and ρ is the density).The class of problems is generally referred to as the Richtmyer-Meshkov (RM) instability [1] [2].One of the eventual goals in investigating RM instability is to shed light on the resultant mixing.Mixing is of contemporary interest in many fields of research, among which are the inertial confinement fusion [3], the fuel mixing in a Scramjet [4], and the explosion of supernovas [5].
In the RM instability researches, an interesting configuration is a shock wave interacting with a cylindrical interface (circular interface in two dimensions) between two fluids.When a planar shock wave impacts on a heavy gas (e.g., SF 6 ) cylinder around by an ambient gas (e.g., air), a shock wave is reflected and a refracted shock wave transmits into the heavy gas cylinder.Because the heavy gas acoustic impedance exceeds that of the ambient gas, the refracted shock is slower than the incident shock wave, and a convergent shock refraction pattern occurs.Because of this and the curvature of the cylinder, the transmitted shock focuses at the downstream vertex.This focusing will induce a pressure rising that eventually leads to a cusp-like protrusion [6].The baroclinically deposited vorticity due to the shock-cylinder interaction stretches and distorts the interface and rolls up into a counter-rotating vortex pair.Then, the evolution of the interface will be dominated by the vortex pair and falls into a self-similar regime, which is in some ways similar to the roll-ups of a vortex sheet [7]- [9].
The propagation of a shock wave in an inhomogeneous medium, and the response of the medium to impulsive acceleration are of fundamental interest in astrophysical systems.The evolution of the interstellar medium in spiral galaxies is significantly influenced by the strong shock waves generated by supernovae explosion [10].Moreover, it is well-known that the generated shock waves remarkably alter the morphology of the cloud (a region of higher density).The bright eastern knot of the Puppis A supernova remnant results in a distorted shock front due to a cloud-shock interaction as seen in images from the Chandra X-ray telescope [11].In this paper, the shock-cylinder interaction is a particularly simple configuration to investigate the problem of shock accelerated inhomogeneous flows.
There are many experimental and numerical researches of the shock-cylinder interactions which concentrate on different subjects.Haas and Sturtevant [12], Jacobs [13] and Tomkins et al. [14] studied this problem experimentally for exploring the wave patterns and the distortion of volume in the RM instability, studying the effects of centrifugal force and viscosity, and investing the mixing mechanisms in a shock-accelerated flow, respectively.Based on the experiments performed by Haas and Sturtevant, Picone and Boris [15] studied the early and the late time phenomena of these experiments numerically.Additionally, Quirk and Karni [16] also simulated the same experiments with the concentration on the early stages of the shock-cylinder interaction.Recently, in light of the experiment performed by Tomkins et al., Weirs et al. [17] investigated the three dimensional effects and Shankar et al. [18] studied the effect of the tracer particle in the shock cylinder interactions.
As be mentioned just, there were some studies on the shock accelerated heavy gas cylinder, however, in the aforementioned papers no attention has been paid on a perhaps existing scaling law in contrast with the single mode or multi mode RM instability [19]- [21].Moreover, most numerical studies were performed for validating codes with comparing to the experimental data.In this work, we systemically and meaningfully analyze the obtained numerical results.Comparing with the results from the prediction of vorticity dynamics, we want to present the behavior of the characterized variables in the evolution of the shock accelerated cylinder.Through studying the evolution of the integral length of the shock-accelerated cylinder, we shed light on the scaling law of the growth of the normalized length follows in the self-similar regime.

Numerical Methods and Initialization
This paper applies our large eddy simulation code MVFT (multi-viscous flow and turbulence) [22] to simulate the multi-viscosity-fluid and turbulence.The code MVFT is used for the compressible large eddy simulation that was developed by Institute of Fluid Physics at the China Academy of Engineering Physics.MVFT can be used to simulate multi-component flows, and compute shocks, contact discontinuities and material interfaces at high accuracy.It splits the flow into an inviscid flow and a viscous flow by using an operator splitting technique, where the former is computed by employing the piecewise parabolic method with a third-order Godunov scheme and the latter is calculated by utilizing a central difference scheme in conjunction with a second-order Runge-Kutta method for the time integration.MVFT applies based on the piecewise parabolic method [23] to interpolate physical quantities, the Vreman [24] subgrid scale eddy viscosity model to conduct large eddy simulation, and to solve Navier-Stokes equations.
The initial conditions are significantly important in numerical simulations, especially for the membrane less RM instability researches.Initially, a sharp interface [25] [26] was used but leaded to an ill-posed simulation compared with the experiments.This is because the interspecies diffusion would not be deniable in the membrane less technique [27] which is used to form the interface between two fluids in experiments.More concentrations should be required on the interfacial diffusion.Consequently, an interfacial transition layer with finite thickness [28] was introduced to characterize the diffusive interface.Recently, a well-characterized initial concentration profile [17] was specified with experimentally measured data.And the contour map of scalar mass fraction and dissipation rate obtained from numerical calculations with the profile were in a good agreement with experimental data qualitatively.For we will simulate the same experiment, the profile is chosen as our initial concentration profile.
In the present simulations, the initial conditions were adapted to the Mach 1.2 shock tube experiment of Tomkins et al. [14] which was performed at Los Alamos National Laboratory (LANL).Figure 1 shows a schematic of the computational domain.To match the shock tube test section dimensions, the computational domain has ( ) The inflow/outflow boundary condition is used on the left and right x boundaries.At the upper y boundary, the reflecting boundary condition is used and the symmetry condition is enforced at the lower y boundary.Three different grid resolutions with x y ∆ =∆ 50 μm, 30 μm and 15 μm = are used for grid convergence study.In the simulations, the initial concentration profile of the cylinder in [17] is used.The simulations are run at a CFL number of 0.2.The main gas parameters are presented in Table 1.

Results
Figure 2 depicts the snapshots of the time evolution of the cylinder, which shows volume fraction maps that are corresponding to the densities for seven times after shock passage on three different grid resolutions.Counter rotating vortex pair formed at the interface owing to the baroclinic vorticity deposited stretches and rolls the interface in ward the vortex cores after the shock collides.At about 220 μs , the vortex pair starts to roll up into the vortex core, and then falls into a self similar state.In Figure 2 one also can see a cusp-like protrusion [6] produced by shock refraction in the simulations with finer grid.This demonstrates that finer grid is needed in the shock-cylinder simulations for better capturing more fine-scale structures.
Figure 3 presents the negative circulation evolution over time of the flow field.Once a straight vortex sheet rolls up, the circulation of each branch is either positive or negative.For the roll-ups of the straight vortex sheet, we only consider the negative circulation here.The amount of circulation deposited during the interaction of the    shock wave with the interface is calculated using a path integration of velocity where V is the velocity vector and s is the path.From Stokes theorem, the circulation is a measure of the vorticity over an area A .In the two dimensional flow, the circulation is ( ) where z ω is the vorticity component which is perpendicular to the x-y plane, u and v are x and y components of the velocity respectively.The vorticity generated by a shock wave propagating through a circular cross-section has been studied by Picon and Boris [29].They gave the magnitude of the vortex strength or circulation, where 2 V is the flow velocity behind the shock in the laboratory frame, W is the shock velocity, 0 R is the ra- dius of the cross-section of the cylinder, ρ ∞ is the ambient density, and c ρ is the density of the gas in the cy- linder.For an initial temperature of 298 K and pressure of 0.8 atm, the resulting 1D gas dynamic velocities are the following, 2 105.6 m s V = , and 414.75 m s W = .In the simulations, the effective radius 0 R is 2.57 mm, the ambient density ρ ∞ is the density of air, 3 0.95 kg m , and c ρ is the density of the cylinder 2.85 kg/m 3 .
From these data, one can get the amount of circulation is 2 0.52 m s − .From Figure 3, one can see that the simulation data is shown a good agreement with the value.
Additionally, one can calculate the vortex strength by applying the approximate model of the compact vortex.When a shock impacts on a heavy gas cylinder, the cylinder will stand here relative to the ambient gas because of its inertia.Then the cylinder will have a velocity 2 V with the opposite shock direction relative to the back- ground.There is a discontinuity of the tangential velocity at both edge of the cylinder, so the cylinder could be considered as a two dimensional vortex sheet.From the vorticity dynamics, one can calculate the magnitude of the circulation of the compact vortex rolled up from the vortex sheet 2 0 2V R Γ = [30].One can obtain a value of 2 0.54 m s − which is also in a good agreement with the numerical results.In two dimensional case, as analogous to center of mass, one can define the coordinate of center of vorticity, The time evolution of the center of vorticity is presented in the Figure 4.One can see that the x-component velocity of center of voriticity is 85.6 m/s.In [15], Picone and Boris gave an equation to compute the perpendicular distance d of the vortex core to the y axis or the half vortex spacing, ( ) Here, cm V is the velocity of center of vorticity in x direction.Then, one can get the value of d, 2.3 mm.From Figure 4, the value is in a good agreement with the center of vorticity in y direction.
The vorticity distribution along y direction computed on the fine mesh resolution at seven different times is shown in Figure 5.The vorticity trends to get together inward and extend outward as time goes by.Because of roll-ups, there are some troughs appeared in the vorticity distribution.After about 400 μs , the vorticity is almost fixed at the outer edge as time elapses, although it also moves inward to the y axis.The behavior is consistent with the evolution of the center of vorticity in y direction which decreases after about 400 μs shown in Figure 4.
Considering the mole fraction ( ) , , X x y t , the spatial maximum of the mole fraction in the stream wise direction gives  .In Figure 6, the log-log plot of normalized length vs. time shows that from about 200 μs to 650 μs , there is a linear relationship.Then, a power law of 0.43 t for the normalized length is obtained.

Discussions and Conclusions
To the authors' knowledge, there hasn't been an investigation of the power law of the length for a shock accelerated gas cylinder, although the power law followed by the mixing zone width has been investigated in the single mode or multi-mode Richtmyer-Meshkov instability.Dimonte and Schneider [19] reported a power law t θ Sohn [21] predicted that the bubble fronts follow a power law with θ in the range of ( ) 0.3 -0.35 0.02 ± by applying a quantitative model of bubble completion.For the roll-ups of vortex sheet, the power law of r t θ with 2 3 θ = has been gained [7], where r is the polar radius of the spiral.Interestingly, it is seen that the power law exponent of the normalized length obtained here is almost same to that of the bubble fronts in [20], although it's quite different from the others.
In conclusion, we studied the evolution of a shock impinging heavy gas cylinder and the growth of the normalized length.A self-similar behavior is shown from the snapshots of the evolution of the gas cylinder.The two dimensional numerical results of the negative circulation, the center of vorticity, and the vortex spacing are in a good agreement with the results obtained from the analyses of Picon and Boris [29].Moreover, the normalized length obeys a power law ~tθ η with 0.43 θ = . Whether the power law is a universal property for a shock accelerated cylinder is under investigation.More investigations of a shock accelerated cylinder with different radii, density ratios, and Mach numbers will be performed in future.

Figure 1 .
Figure 1.Schematic diagram of the computational domain.

Figure 2 .
Figure 2. Image sequence of SF 6 volume fraction at t = 130, 220, 310, 400, 560, and 650 μs after shock impingement on the coarse (left column), medium (middle column), and fine (right column) mesh resolutions.The shock has traversed the cylinder from top to bottom.

Figure 3 .
Figure 3. Negative circulation as a function of time on three different grid resolutions.The short dashed, short dotted and solid lines correspond to the results computed on the coarse, medium, and fine grid resolutions respectively.The dashed line is the theoretical prediction for the negative circulation in [29].

Figure 4 .
Figure 4. Center of vorticity in x direction (solid line) and y direction (shot dot) as a function of time computed on the fine grid resolution.The dashed line is obtained from linear fitting.The dash dot line corresponds to the theoretical prediction for the half of vortex spacing in [29].

Figure 5 .
Figure 5. Vorticity distribution along y direction at seven different times on the resolution of 15 μm x y ∆ =∆ = .

Figure 6 .
Figure 6.Log-plot of normalized length vs. time computed on the coarse (shot dashed line), medium (short dotted line) and fine (solid line) mesh resolutions.The dashed line corresponds to the power law of 0.43t .The scattered data correspond to the experimental data obtained from[14].

Table 1 .
Properties of air and SF 6 gases.