# Shock tube

### NUMERICAL MODELING OF A REFLECTED SHOCK WAVE DURING SHOCK-TUBE IGNITION DELAY TIME EXPERIMENTS AT PRACTICAL CONDITIONS

### Abstract

A numerical model for reflected shock wave during ignition shock-tube ignition delay time experiments is described, in which all the processes of bifurcation, wall effects, heat transfer, and shock flame interactions are incorporated. The model is being developed for both cold flow and reactive flow using methane and hydrogen mixtures. The numerical method is first applied to the problem of one-dimensional, no expansion with slip condition for a normal shock wave formation in a closed shock tube then bifurcation effects are incorporated. A three model is being developed and compared to an axisymmetric model. The effect of wall disturbance and heat transfer with and without boundary layer effects are being investigated. The simulation is being performed for both cold and reactive flow. The model is valuable to the interpretation of the experimental results and to the understanding of various phenomena occurring in shock tube experiments. The model should be able to predict temperature, pressure, and density profiles behind the reflected shock wave. The objective of this study is to perform numerical simulations of the multidimensional physical processes occurring during shock tube experiments. First, a general overview on shock tube experiments is provided. Then a combustion time model developed from ignition delay time data is introduced. Last, an approach on the basic concepts necessary to setup the numerical model to simulate the different phenomena induced by ignition behind reflected shock wave is outlined. The multidimensional, time-dependent numerical simulation of unsteady combustion should be able to resolve all of the relevant scales, ranging from the size of the system to the reaction zone scale.

### GENERAL INTRODUCTION

The shock tube has found widespread use as an experimental device in which to investigate chemical kinetic behavior in reactive gas mixtures. A shock tube consists of a high-pressure driver section and low-pressure driven section initially separated by a diaphragm. The driver section is pressurized to energy high enough to cause the diaphragm to rupture and as a result, a shock wave is generated and travels down the driven tube. Simultaneously, an expansion fan propagates through the high-pressure side. Both waves reflect off the shock tube endwalls. Of interest herein is the endwall region behind the reflected shock wave where kinetic experiments take place. Ideally the flow properties behind the reflected shock are uniform and do not vary with time. In real shock tubes however, there is a slight deviation from this assumption due to viscous and non-ideal effects. As a result, some non-ideal effects occur that impact the test conditions, the magnitude of which depends on the tube diameter, pressure, and shock Mach number.

### The Shock tube Facility

A shock tube consists of a rigid cylinder in which a gas at high pressure, called the driver gas, is initially separated from a gas at lower pressure, called the test gas, by a diaphragm. The high and low pressure regions can alos be referred to as compression and expansion pressure chambers respectively (Wright, 1961). When the diaphragm is suddenly burst, an incident shock wave propagates through the test gas raising the temperature and pressure. The incident shock arrives at the end wall and reflects back raising the temperature even further to start energy release and therefore chemical reaction commences which leads to ignition. As the shock wave moves through the test gas, a rarefaction wave moves back into the high-pressure gas at the speed of sound. The test gas and the driver gas make contact at the “contact surface”, which moves along the tube behind the shock front. The end of the experiment is dictated by the intersection of the reflected shock and the expansion fan. Typical test times are on the order of few milliseconds. The shock wave will be unattenuated with distance, provided the tube is of constant area. This property allows for a controlled shock wave, which makes it an invaluable tool in many investigations (Wright, 1961). 4 shows the ideal movement of the shock front, the contact surface, the rarefaction wave and the reflected shock wave in a distance-time diagram. The unperturbed, low-pressure test gas is given by by the subscript 1, and the initial temperature and pressure in this region are denoted as p1 and T1, respectively. The region between the shock front and the contact surface is denoted by 2; the region between the contact surface and the rarefaction wave by 3. The initial conditions on the high-pressure side are given the subscript 4. If the shock wave is allowed to undergo reflection at the end of the tube, the pressure conditions in this region are given the subscript 5.

There exist two shock tube facilities in the research laboratory of Dr. Eric Petersen at Texas A&M University. Ignition delay time data can be generated in both test facilities. The smaller shock tube facility generates pressures up to 10 atm, while the bigger facility generates pressures up to 100 atm. The 100 atm shock tube has a driver length of 3.5 m, and a driven length of 10.7 m. The inner diameter is 16.2 cm. The facility is described in details by Petersen et al (Petersen E. L., 2005). The smaller shock tube is made of stainless steel and has a total length of 6.1 m. The square, helium-driven section is 4.3 m long and has a 10.8-cm height. The circular driver tube is 1.8 meters in length with a 7.60-cm inner diameter. A schematic of the shock tube is given in Fig. 1.

Target reflected-shock temperatures and pressures are achieved through the use of 1-mm Lexan diaphragms. Ignition measurements are performed in the reflected-shock region. Incident-shock velocities are measured at four different axial locations along the driven tube by using four PCB 113A pressure transducers in conjunction with four 120-MHz Fluke model PM6666 time-interval counters. Temperatures and pressures in the reflected-shock region are determined from the standard 1-D shock-tube relations and the Sandia thermodynamic database. Ignition delay times are measured from the endwall by monitoring the endwall pressure signal through a PCB 113A pressure transducer. The endwall emission trace is also used as a guide for ignition delay time determination and is acquired by monitoring the CH* chemiluminescence through a 430±5 nm narrow band filter with a Hamamatsu 1P21 Photomultiplier tube (PMT) in a custom-built housing. Optical ports located on the side of the shock tube allow for additional non-intrusive optical access. 2 shows sample endwall pressure and endwall emission traces. The sidewall pressure and sidewall emission traces were also recorded but were not used to determine the ignition times (Fig. 3). The data acquisition system comprises of two 16-bit 10-MHz computer oscilloscope boards with a total of four channels.

### Shock-Tube Non-Idealities

There have been a number of experimental studies investigating the mechanisms responsible for non-ideal behavior in shock-tube facilities. These mechanisms include the non-ideal rupture of the diaphragm,1-6 (1) (2) (3) (4) (5) (6)reflected shock/boundary layer interactions,7-14 (7) (8) (9) (10) (11) (12) (13) (14)driver gas contamination,15,16 (15) (16)contact surface instabilities,17-22 (17) (18) (19) (20) (21) (22)and thermal boundary layer effects.23-25 (23) (24) (25)

Also, theoretical and analytical models have been developed that provide a reliable understanding of the non-idealities in the shock tube. In particular, Mark was the first to establish a theoretical treatment of the reflected shock/boundary layer interaction and the mechanism responsible for wall jetting of driver gas through the bifurcated structure.7 (7) Since then, several empirical models have been developed in an attempt to quantify the disturbance level of the bifurcated zone.26,10,12 (26) (10) (12)Davies and Wilson27 (Davies, Influence of Reflected Shock and Boundary-Layer Interaction on Shock Tube Flows, 1969)and Stalker and Crane28 (28)used Mark's theory and developed analytical models for predicting the premature arrival of driver gas to the endwall region. Numerical simulations have also provided reliable information about the contamination process which is of great concern in high-enthalpy shock tunnels in particular.27,28 (Davies, Influence of Reflected Shock and Boundary-Layer Interaction on Shock Tube Flows, 1969) (28)

Ideally, the endwall region behind the reflected shock wave is assumed to be isothermal before chemical reaction takes place. However, in real shock-tube experiments and especially when the test times are relatively long, the isothermal assumption becomes invalid, and heat losses from the hot gas and the cold wall become important and should be accounted for in shock-tube modeling studies. Furthermore, the effect of the thermal boundary layer on the reflected-shock structure and the flow properties near the endwall region have been analyzed theoretically29-32 (29) (30) (31) (32) and numerically based on the gas kinetic models (BGK)33 (33)and the 1-D Navier-Stokes equations.34,35 (34) (35)

Unfortunately, the shock tube is a transient test facility with unsteady and highly nonlinear physical processes that cannot be accurately modeled with simple 1-D models, but the investigation of complex phenomena associated with the shock tube can be enhanced considerably when done in concert with multi-dimensional numerical simulations. Such simulations are made possible with Computational fluid Dynamics (CFD) codes and the increased computational resources which have become quite powerful tools to highlight the flow physics in multi-dimensional complex and transient flow fields.

To that end, shock-tube non-idealities have been investigated during the last several decades. Multi-dimensional simulations have been performed to model the non-ideal rupture of the diaphragm to quantify the shock speeds, the structure of the developing flow downstream the diaphragm and the contact surface shape. For example, the opening of the diaphragm has been modeled as a slit in two dimensional simulations,6,36 (6) (36)versus an iris in axi-symmetric simulations.37-40 (37) (38) (39) (40)The reflected-shock bifurcation phenomenon has been simulated in a number of numerical studies14,41-51 (14) (41) (42) (43) (44) (45) (46) (47) (48) (49) (50) (51)and also has served as a test case for numerical method validations.51-56 (51) (52) (53) (54) (55) (56)Multi-dimensional simulations have also aimed at modeling the interaction between the reflected shock and the contact surface in an attempt to predict the driver gas contamination in high-enthalpy shock tunnels in particular.38-40,46,57-61 (38) (39) (40) (46) (57) (58) (59) (60) (61)

### Motivation

While the studies mentioned above have given insight into the flow evolution in shock tubes and the various non-ideal phenomena that affect the test time and conditions, the inability to model the whole test facility geometry due to computational limitations necessitates various assumptions and consequently, the majority of the previously reported simulations focused on certain parts of the shock-tube facility such as the endwall region with the upstream inflow conditions derived from the Rankine-Hugoniot relations and the initial conditions set just before shock reflection upstream of the endwall. When viscous effects are taken into consideration, the initial flow field behind the incident shock is usually estimated from the boundary layer theory of Mirels.62,63 (62) (63)Therefore the results depend greatly on the assumptions associated with modeling only part of the facility. Of particular interest in this study is the development of a fluid mechanics model of the reflected-shock process that can be applied to experimental conditions routinely seen in high-pressure chemical kinetic and ignition delay time experiments in undiluted fuel-air mixtures. With such a model, parametric studies can be performed that couple the shock-tube fluid mechanics with the chemical kinetics in an attempt to characterize the extent of non-ideal behavior in such experiments.

In this respect, accurate simulations of the complex spatial and temporal inter relations in the shock-tube flow are more likely to be achieved through modeling of the complete shock-tube geometry rather than just the endwall region of the shock tube. This approach is a challenging task because it requires the use of large mesh sizes in addition to time marching the solution over a large number of small time steps that allow the resolution of the flow-physics. Such large-scale computations can be made possible through the use of parallel processing.

### Study Outline

As a first step towards modeling the complex mechanisms responsible for the non-uniform conditions and the reduced test times in the shock tube, the first half of this study assumes the flow is inviscid and focuses on developing a robust and accurate model capable of reproducing the major flow phenomena in the shock tube from the propagation of the shock wave and contact surface to the reflection of the expansion fan and incident shock. The inviscid model should serve as a baseline for the rather more complex viscous model which necessitates increased computational efforts. The modeling of the complex flow structure in the shock tube by solving the unsteady Euler equations under the assumption of inviscid flow have been successfully investigated in 1-D,64-66 (64) (65) (66)2-D,6,36,67-69, (6) (36) (67) (68) (69)axi-symmetric,37 (37) and 3-D70,71 (70) (71)simulations. In addition, the axi-symmetric approach has proven to be an appropriate representation of the cylindrical geometry of the shock tube.37-40,46,49,57-61 (37) (38) (39) (40) (46) (57) (58) (59) (60) (61) (49)

However, when shock-tube non-idealities need to be modeled to quantify their impact on the flow uniformity in the test region, viscous effects have to be taken into consideration. In this respect, the second half of the study introduces a viscous model of the shock tube which simulates the major non-ideal mechanisms responsible for non-ideal behavior, including the reflected shock/boundary layer interaction or bifurcation, reflected shock/contact surface interactions, and heat transfer from the hot test gas to the cold shock-tube walls.

As a first step towards modeling the non-idealities in the shock tube, the boundary layer is assumed as laminar, and a conjugate heat transfer model is incorporated to account for the heat losses from the test region hot gas to the cold shock-tube walls. The accuracy of the simulations is enhanced by modeling the boundary layer in the shock tube as turbulent and the k-e realizable turbulence model was used. The turbulent boundary layer mesh was carefully treated, and the enhanced wall treatments were utilized to ensure a y+ value close to 1 and to guarantee the flow structure is resolved all the way to the viscous sublayer.

The results from the axi-symmetric simulations of the shock-tube model are presented in this paper. The simulations were carried out with the commercial CFD solver FLUENT using the coupled explicit density solver and the finite-volume approaches. The accuracy, efficiency, and stability of the numerical model are investigated with 1st order, 2nd order, 1st/2nd order blending, and 3rd order MUSCL schemes. Adaption was used to assure mesh refinement in high-gradient regions to accurately resolve the shock and contact discontinuities. The shock-tube geometry is represented by a structured axi-symmetric mesh.

A general overview of the model setup and governing equations is first presented. Then descriptions of the model geometry, mesh, boundary and initial conditions are presented. Provided next is the simulated inviscid solution of the temporal and spatial variation of the flow properties which are validated with the 1-D inviscid theory and with experimental data. Then the viscous solution is presented and the modeling and the bifurcation phenomenon due to shock wave-boundary layer interaction is simulated. Finally, a discussion of the results and the conclusions are given.

### COMPUTATIONAL MODEL

The simulations were carried out with the commercial CFD solver FLUENT. The model equations are discretized in space and time following the control volume approach and utilizing a density-based explicit solver. The flow domain is represented with a structured mesh of hexahedral cells. A grid adaption tool is used to resolve regions with the steepest gradients. The AUSM+ flux vector splitting scheme is used to compute the flux vectors. The convective terms are discretized in space following three schemes, 1/2nd order blending upwind, 2nd order upwind, and 3rd order MUSL schemes to investigate the impact of increasing scheme resolution on the accuracy and stability of the solution. An explicit time-stepping integration was performed using a four-stage Runge-Kutta scheme for unsteady flows. The time step was set by the Courant-Friedrichs-Lewy stability limit between 0.8 and 1.

### Grid

The computational domain represents the entire geometry of the high-pressure shock tube test facility at Texas A&M University described in detail in de Vries et al.72 (72) The shock tube consists of a 2.46-m long driver section with an internal diameter of 7.62 cm and a 4.72-m long driven section with an internal diameter of 15.24 cm. The shock-tube side wall and endwall thicknesses are 1.27 cm and 2.54 cm, respectively. Due to the axial symmetry of the shock tube, the flow domain is modeled with an axi-symmetric structured mesh as shown in Fig. 1. The axi-symmetric approach is appropriate for the cylindrical geometry of the shock tube and is sufficient to render an accurate description of the real flow configuration. The boundary conditions of the shock-tube model when the conjugate heat transfer model is turned on are given in Fig. 2.

### Finite Volume Method

The density-based coupled-explicit algorithm is adopted with double precision. This algorithm solves the governing equations of continuity, momentum, energy and species transport simultaneously as a set of equations. The finite-volume-based discretization approach is used to solve numerically the Euler equations. This is accomplished by first dividing the domain into discrete control volumes, then integrating the governing equations on the individual control volumes to construct algebraic equations for the discrete dependent variables and last linerazing the governing equation to produce a system of equations for the dependent variable. The governing equations can be illustrated by considering the unsteady conservation equation for transport of the scalar quantityφ. This is demonstrated by the following equation written in integral form for an arbitrary control volume V as follows,

∂ρφ∂tdV+ρφv.dA=Γφ∇φ.dA+SφdV

(1)

Where, ρ is the density, v is the velocity vector, A is the surface area vector, Γφ is the diffusion coefficient for φ, ∇φ is the gradient of φ, and Sφ is the source of per unit volume. Discretization of Eq. (13)on a given cell yields

∂ρφ∂tV+fNfacesρfvfφf.Af=fNfacesΓφ∇φf.Af+SφV

(2)

Where Nfaces are the number of faces enclosing the cell, φf is the value of φ convected through face f, ρfvf.Af is the mass flux through the face f, Af is the area of face f, ∇φf is the gradient of φ at face f, and V is the cell volume. Linearization of equation is performed as such,

apφ=nbanbφnb+b

(3)

Where the subscript nb refers to neighbor cells, and ap and anb are the linearized coefficients for φ and φnb.

### Spatial Discretization

The AUSM+ flux vector splitting scheme was used to compute the flux vectors. AUSM stands for Advection Upstream Splitting Method first introduced by Liou and Steffen89, accurately captures shock and contact discontinuities by providing an exact resolution, free of “Carbuncle” phenomena and oscillations at stationary and moving shocks, and preserves the positivity of scalar quantities. The AUSM scheme was substantially improved to yield the generalized Mach number-based convection and pressure splitting functions. The new scheme termed AUSM+ is more robust and avoids using an explicit artificial dissipation.86,90

### Order of Accuracy

First order can yield better convergence than higher order schemes especially in complex flows where shocks and discontinuities are present, however, especially when turbulence is dominant, first-order discretization is not recommended due to the increased numerical discretization error introduced by this scheme. First-order accuracy is obtained by simply setting the face value φf equal to the cell-center value of φ in the upstream cell as shown in Eq. (16). fou stands for first order upwinding.

φf,fOU=φ

(4)

In general, when shock discontinuities are present in the flow, it is almost impossible to achieve a stable solution, free of unphysical numerical oscillations and nonlinear instabilities without introducing some numerical dissipation. However, numerical dissipation can be reduced by increasing the mesh resolution at the shock and contact discontinuities by applying the adaptive grid refinement approach.

High resolution schemes such as 2nd order or higher schemes can yield more accurate results than the first order scheme. If second-order accuracy is desired, quantities at cell faces are computed using a multidimensional linear reconstruction approachby applying Taylor series expansion about the cell centroid. The face value φf is computed using the following equation:

φf,SOU=φ+∇φ.r

(5)

Where SOU stands for second order upwinding, φ is the cell-centered value and ∇φ is its gradient in the upstream cell, r is the displacement vector from the upstream cell centroid to the face centroid. For second order schemes, oscillations creep up, notably at the discontinuities. A Total Variation Diminishing concept is applied and the gradient ∇φ is limited so that no new maxima or minima are introduced.

Another effective method to avoid the oscillations in the solution is by adding some diffusion to the 2nd order scheme by blending the second order φf,SOU and first order φf,fOU approximations as such,

φblending=φf,fOU+β(φf,SOU-φf,fOU)

(6)

Where β is a blending factor whose values lie between zero and one, β = 0 reduces the scheme to 1st order upwind, and β = 1 brings it back to pure 2nd order upwinding providing satisfactory results. Usually small amounts of the first order scheme (10-20%) is sufficient to get rid of the oscillations, and the accuracy is nearly as good as with 2nd order accurate scheme and the stability is nearly as good as the 1st order scheme.

Third order MUSCL Scheme79 (Monotone Upstream-Centered Schemes for Conservation Laws)78 is achieved by blending a central differencing scheme and second order upwind scheme as shown in Eq. (19).

φf,MUSCL=θφf,CD+(1-θ)φf,SOU

(7)

Where φf,CD is determined using central differencing and φf,SOU is computed using the second-order upwind scheme as described in Eq. (17). Spatial accuracy is improved with the MUSCL scheme by reducing numerical diffusion, however the implemented scheme in the CFD solver FLUENT does not contain any flux-limiter and therefore produces undershoot and overshoot in the solution.

Temporal Discretization

Time marching is performed by evaluating the scalar φ at the current time level

φn+1-φn∆t=F(φn)

(8)

And φn+1 can be explicitly expressed in terms of the existing solution values φn as such,

φn+1=φn+∆t F(φn)

(9)

Here, the time step ∆t is restricted to the stability limit set by the Courant-Friedrich-Lewy condition CFL as shown in Eq. (22). To maintain time accuracy of the solution, the explicit time stepping employs the same time step in each cell of the domain. The use of explicit time stepping is fairly restrictive. It is used primarily to capture the transient behavior of moving waves, such as shocks, because it is more accurate and less expensive than the implicit time stepping methods in such cases. A four-stage Runge-Kutta scheme for unsteady flows is used by the CFD solver. The stability limit that was set between 0.8 and 1. The time step ∆t is given by,

∆t=2 CFL.VfλfmaxAf

(10)

Where V is the cell volume, Af is the face cell, and λfmax is the maximum of the local eigenvalues.

### Dynamic Grid Adaption

Adaptive grid feature was used to cluster the grid around the shock and contact discontinuities. As with any numerical approximation of solution discontinuities, shocks will be somewhat smeared by grid resolution and numerical diffusion, adaption is an effective tool that can be used to resolve the discontinuities and reduce the numerical error in the digital solution with minimal numerical cost.91,92 The adaption feature can also be used to achieve a grid-independent solution without regenerating the mesh.

Adaption parameters can be tuned to effectively capture shock and contact discontinuities. Coarsen and refine thresholds are adjusted to achieve the desired level of adaption.

Adaption is performed by refining the mesh near high gradient regions and coarsening the mesh wherever else needed, thus providing better resolution of shock and contact discontinuities. The refine and coarsen limits are determined based on the normalized values of the density gradients. In this approach, the Euclidean norm of the gradient of the selected solution variable is multiplied by a characteristic length scaleand the gradient function has the following form,

ei1=Acellr2∇f

(11)

Where ei1 is the error indicator, Acell is the cell area, r is the gradient volume weight with r = 1 corresponding to full volume weighting, and ∇f is the Euclidean norm of the gradient of the desired field variable f. The normalized values are obtained by scaling the values of ei1 by their maximum value in the domain as is shown in Eq. (24) so that the adaption rang, refine and coarsen thresholds, is between [0, 1]. The refine and coarsen threshold values depend on the strength of the shock and contact discontinuities.

ei1maxei1

(12)

Dynamic grid adaption based on gradients of both density and pressure was performed every five iterations.

### INVISCID SOLUTION

### Introduction

As a first step towards modeling the complex mechanisms responsible for the non-uniform conditions and the reduced test times in the shock tube, the first half of this study assumes the flow is inviscid and focuses on developing a robust and accurate model capable of reproducing the major flow phenomena in the shock tube from the propagation of the shock wave and contact surface to the reflection of the expansion fan and incident shock. The inviscid model should serve as a baseline for the rather more complex viscous model which necessitates increased computational efforts. The modeling of the complex flow structure in the shock tube by solving the unsteady Euler equations under the assumption of inviscid flow have been successfully investigated in 1-D,64-66 (64) (65) (66)2-D,6,36,67-69, (6) (36) (67) (68) (69)axi-symmetric,37 (37) and 3-D70,71 (70) (71)simulations. In addition, the axi-symmetric approach has proven to be an appropriate representation of the cylindrical geometry of the shock tube.37-40,46,49,57-61 (37) (38) (39) (40) (46) (57) (58) (59) (60) (61) (49)

### CONCLUSIONS

In order to understand the phenomena occurring during shock-tube ignition experiments, a numerical developed is being developed in order to simulate the different shock-physical events interactions which are not obviously detected in laboratory experiments. A numerical model will not only provide a deep insight on the physical and chemical processes taking place in a shock tube, but will also serve as a tool to simulate conditions and parameters that are not possible to determine by the mean of experiments alone. This project will focus on addressing the issues provided in the outline section. The modeling of the reflected shock with the different fluid and heat dynamics will be first based on basic one-dimensional principles all the way to three-dimensional simulation with added reaction. The large disparity in the spatial and temporal scales that can vary by up to 6 orders of magnitude will further add to the complexity of the modeling. Another complexity factor is the incorporation of the detailed chemical reaction mechanism to the numerical model. However, the chemistry model proposed herein was developed from the detailed mechanism and is modeled as single step time model. This model, when incorporated with the Eddy-Breakup model produces results with accuracy and in the most computationally efficient way. As part of setting up the numerical model, boundary conditions will be identified and grid models will be generated. The numerical code is based on a CFD package that was developed by Flow Parametrics. The simulations results achieved by the Oran group will be used as a validation tool for this study.

### Recommendations

APPENDIX A:

APPENDIX B:

LIST OF REFERENCES