1D study of the detonation phenomenon and its influence on the interior ballistics of the combustion light gas gun

Fei Zhou , Ning Liu ,*, Xing-yn Zhng , Xiu-dong You

a School of Mechanical Engineering, Nanjing University of Science and Technology, Nanjing 210094, Jiangsu, China

b Troops 91388 of PLA, Zhanjiang 524022, Guangdong, China

Keywords:Riemann solver WAF Combustion light gas gun Detonation

ABSTRACT One-dimensional simulations with a detailed hydrogen-oxygen reaction mechanism have been performed to investigate detonation phenomenon in a combustion light gas gun (CLGG). Convection fluxes of the Navier-Stokes equations are solved using the WAF(weighted average flux)scheme HLLC Riemann solver.A high resolution fifth-order WENO scheme for the variable extrapolation at the volume interface and a fourth-order Runge-Kutta scheme for the time advancement are used. Validation tests of the stoichiometric hydrogen-oxygen deflagration to detonation transition process shows good agreement between the computed results and the analytical and documented solutions, demonstrating the reliability on the detonation simulation of the current scheme. Simulation results of the interior ballistic process of the CLGG show that the flame propagation experiences three distinct stages. The blast detonation wave causes a high-pressure shock and hazardous oscillations in the chamber and makes the projectile accelerate with fluctuations, but has only a small improvement to the muzzle velocity.

1. Introduction

The combustion light gas gun(CLGG)is an advanced launching system with potentially extremely high muzzle velocity, in which the projectile is propelled by the combustion expansion of a low molecular weight combustible gas mixture[1].The combustible gas mixture is oxygen and hydrogen. A high charge density of the hydrogen-oxygen propelling charge makes the flame propagation accelerate fast and may trigger a deflagration to detonation transition (DDT).

Numerical research attempts to clarify the DDT mechanism have a long history that begins with Zel'dovich [2], where a onedimensional single reaction step model was used. Their simulation results reveal that detonation can be generated when a flame propagates through an unburnt gas mixture with a temperature gradient.Oran and Gamezo[3]simulated the DDT process induced by a localized hot-spot in the sensitized region; firstly the deflagration flame accelerates into a local explosion through the gradient mechanism, then the eventual detonation is triggered as the local explosion continues accelerating. They pointed out that shockwaves strengthen the turbulent fluctuation of the unburnt nearby flame front, which in turn enhances the thermal exchange and the species transport. Direct numerical simulation (DNS) results of Ivanov et al. [4] proposed a different DDT mechanism in a smooth tube. It was found that a shock wave may propagate overlapped with the reacting region at the end of the fast deflagration stage,as a result,the continuous enhancement between the shock wave and reaction heat release strengthens the flame velocity and the shock pressure exponentially,which triggers the final detonation. Other mechanisms leading to the onset of detonation can also be seen reported [5,6], such as shock focusing and shock induced auto ignition, these cases are usually associated with boundary conditions.

It is clear that 2-D (2-dimensional) and 3-D (3-dimensional)simulations have obvious advantages on the numerical research of DDT phenomenon, due to the expression of important multidimensional details,especially small-scale characteristics around the thin reacting zone. However, the computing cost of a multidimensional simulation is also expensive,which is mainly caused by two reasons. The first one is the size of the computing matrix,which is dominated by the number of the grids. As the grid independency research of Ivanov [7] and Heidari [8] indicates, DDT simulation requires at least 6-8 grids per flame width for acceptable solution convergence or 16 grids for negligibly small numerical error. The grid size for stoichiometric H2-O2combustion simulation should be smaller than Δ=0.0375 mm (about 1/8 of the laminar flame thicknessLf≈0.3 mm)atP0=105Pa,T0=298 K,and a smaller size is necessary as the initial pressure grows. Consequently, millions of grids are usually required for the DDT simulation of medium scale problems. The second one is the number of species transport equations, which is controlled by the reaction mechanism.Research of Liberman et al.[9,10]showed that a singlestep reaction model gives error prediction on the minimal length of the temperature gradient capable to trigger a detonation, hence a detailed or simplified multi-step reaction mechanism with key chain branching reactions is necessary.

To reduce the computing cost, various efforts on the model simplification can be seen reported, such as adaptive mesh refinement method[11]for grid-number reduction,in situ adaptive tabulation [12] for reaction iteration acceleration, artificial thickened flame [13] and LES model [14,15] for coarse grid application.However,as the flame surface becomes thinner at higher pressure,the minimum grid size required decrease, the necessary grid number to solve the CLGG ballistics by a multi-dimensional model based on micron scale grids is tens of millions of magnitudes.Therefore, a reactive 1-D model is more attractive for the detonation influence analysis of the CLGG ballistics.

To meet the shock-capture precision requirements on the simulation of large gradient flow problems, such as dam-break flow, supersonic combustion and DDT phenomenon, various schemes have been proposed, such as Roe with momentum interpolation[16],AUSM+,TVD-MUSCL,WENO,etc.The HLLC Riemann solver is an improved method of the HLL solver by making better characteristic wave assumption,and has been proved to have good adaptability on internal ballistics simulation [17]. The weighted average flux (WAF) approximation method is a second-order extension of the first-order Godunov method, which was first proposed by Toro [18]. When used in conjunction with the HLLC Riemann solver [19], it guarantees a second order of accuracy in time and space with no need for data reconstruction, which has been proved to have good shock-capture accuracy for different problems [20,21].

This work devotes to applying the WAF scheme HLLC Riemann solver to the numerical study of the DDT problem in CLGG.Combining with a detailed 19-step reaction model,a high temporal and spatial resolution 1-D model is performed, DDT phenomenon of high pressure hydrogen-oxygen flame under different initial component ratios is simulated, the influence on the interior ballistics is analyzed.

2. Model description

2.1. Governing equations

For transient compressible reactive flows, the 1-D N-S governing equations with a non-equilibrium chemical reaction mechanism can be written as follows:

In the above equations, the mass, momentum, energy and species transport scalars ρ, ρu, ρEand ρYkare solved simultaneously, while other scalars are solved afterward. Symbol A represents the cross-sectional area of the chamber.

The pressure of the gas mixture is calculated using the Peng-Robinson equation of state:

where the R is the universal gas constant,the formation ofVm，bm，am（T）are taken from the NIST database[22]with molar-weighted average for the mixture parameter calculation.

The Soret form of the thermal diffusion coefficientDT，kcan be written as follows:

whereMw，k,XkandYkare the mole weight,mole fraction and mass fraction of componentk.

The viscosity μ and the specific heatCpof the gas mixture are the molar-weighted average value of each component.The viscosity of componentkis a function of temperature and can be calculated by the Sutherland law with three coefficients,while the thermal conductivity is estimated using the Eucken formula,κT= μk（Cpk+ 0.8R）.

The specific heatCpkand the enthalpy of componentkare calculated by the thermodynamic polynomials taken from the NIST database[22]:

where ΔHkis the enthalpy of formation of componentk.

2.2. Finite rate reaction model

The reaction part of the source termis calculated based on the Arrhenius equation and ΔHkis the enthalpy of formation of componentk, the total generation rate of componentkis defined as:

whereNrandNare the number of the reaction equations and the species,andare the product and reactant stoichiometric coefficients for specieskin reactionr,andare the product and reactant rate exponent for species k in reactionr,Cjis the molar concentration of speciesjand Γrrepresents the effect of third bodies.

The forward reaction rate in Arrhenius form iskf=ATβe-E/（RT）,whereA, β andEare the pre-exponential factor, temperature exponent and activation energy.

The present chemistry model includes 9 species and 19 elementary reaction steps,which is taken from Burke et al.[23]and has good applicability in a wide range of initial temperature,pressure (0.1～≥100 atm) and equivalence ratio (0.1-5) [24,25].

2.3. The WAF method

First, let us consider the 1-D Euler equation by dropping the diffusion and source terms in equation (1), which becomes:

According to the arrangement of the three-wave structure HLLC solver in Toro[18],a non-oscillatory TVD version of the WAF flux at the right interface of cellican be written as:

The Roe average is used for calculating the interface coefficient in the diffusion terms,and then the diffusion terms are discretized by second order central difference scheme. This method is simple and has relatively low precision,however,considering the diffusion terms are usually two or three orders of magnitude smaller than the convection terms,the numerical error is acceptable.

2.4. Time marching scheme

A four-stage fourth-order explicit Runge-Kutta scheme is adopted for time marching, which is given by:

Where α = ［0，1/2，1/2，1］ and β = ［1/6，1/3，1/3，1/6］.

The governing equations are solved using a parallel MATLAB code on an Intel E5-2680V2 workstation.

3. Validation test of the detonation phenomenon of the stoichiometric H2-O2 flame

Simulation of the stoichiometric H2-O2flame propagation with initial temperature and pressure ofT0=298 K andP0=105Pa is performed to validate the model applicability on DDT problems,and the results are compared with the experimental data of Kuznetsov et al. [27].

The premixed mixture is ignited by a small pocket of 2000 K,5 mm fully burnt gas on the left end of the tube,and the length of the computational domain isL=4 m,which is long enough for the flame to accelerate to detonation.To reduce the computation cost,the domain is divided into 16 uneven zones and grid refinement is applied to the regions with flame propagating through.

Grid independency tests of the current model with grid sizes△ =20 μm,10 μm,6 μm,4 μm,3 μm,2.5 μm have been carried out,the predicted flame positions are compared in Fig. 1. Detonation transition occurs atx=925.1 mm,815.6 mm,781.7 mm,744.3 mm,748.4 mm, 747.0 mm respectively. When the grid resolution is coarser than 10 μm, which is about 1/5 flame thickness, the flame acceleration becomes significantly slower and the detonation transition needs a longer distance. The solution tends to converge on 12 and more grids per flame thickness,and the results obtained with the grid resolutions of 3 μm and 2.5 μm are almost the same.

After ignition,the flame propagation undergoes an exponential acceleration and then moderates to a linear acceleration until the flame velocity rises above 700 m/s, as shown in Fig. 2. Typical phenomenon of the formation process of the preheat zone and the precursor shock wave can be observed from the temperature and pressure profiles in Fig. 3, which is caused by the compression heating and the superimposition of the strengthening pressure waves propagating from the flame front. The temperature in the preheat zone before detonation transition is lower than 550 K.

Fig.1. Comparison of the flame front position-time evolution for different grid sizes,△ =20 μm,10 μm, 6 μm, 4 μm, 3 μm, 2.5 μm, which are approximately 2/5,1/5,1/8,1/12,1/16 and 1/20 of the laminar flame thickness, the inset figure is the flame location deviation against △ =2.5 μm.

Fig. 2. The velocity-time dependence of the combustion wave, the fitted curve in the enlargement representing the exponential acceleration stage is approximately 40.337 exp(4.137 t/ms).

Detonation transition occurs attDDT=1.976 ms andxDDT=747.9 mm, the deflagration flame steepens into an overdriven detonation with the maximum pressure of 4.76 MPa and subsequently decays towards a Chapman-Jouguet detonation(DC-J)with computed pressure and velocity of 2.03 MPa and 2771 m/s.The simulation errors on the C-J velocity and C-J pressure against the experimental values in Ref. [27] are -2.4% and 6.8%, while the erro on the DDT location is 6.3%,hence the model proposed in the paper is acceptable for the simulation of detonation problems.

4. Detonation phenomenon in the interior ballistic process of CLGG

4.1. Physical settings and basic assumptions

Fig.3. Temperature(dashed lines)and pressure(solid lines)distribution representing the flame front and the pressure peak formation. Time instants are from t1=0.5 ms,Δt=0.2 ms, tend=1.9 ms.

Fig. 4. Schematic diagram of 45 mm CLGG with 5 ignition points.

Table 1 Initial parameters of a 45 mm CLGG.

Fig. 5. Flame velocity evolution of 3-point ignition and MH2:MO2=6.

Fig. 6. Time-varying temperature distribution reflecting the flame acceleration process.

A schematic view of the 45 mm CLGG is shown in Fig.4,and the initial parameters are listed in Table 1. The chamber is sealed by a gasket ring which breaks above the starting pressure. Interior ballistics with different ignition points and hydrogen-oxygen molar ratios are simulated separately, representing the detonation occasion and the normal combustion process which has no detonation.

In current model,the cross-section areaAis a linear function of the axial location,and the minimum length of the grids is 0.15 μm,which is about 1/16 of the flame thickness. Besides, it is assumed that no mass leaks during the ballistic cycle and the ignition is simplified by patching a small pocket of 2000 K, 4 mm width fully burnt region to the ignition location, which provides a mild flame initiation.

4.2. Detonation phenomenon of 3-point ignition and hydrogenoxygen molar ratio of 6:1

Fig. 5 and Fig. 6 shows the temperature distribution and flame velocity history till fast deflagration formation. The early stage of the flame acceleration is similar to the normal pressure situations.Flame acceleration experiences a period of exponential rising after ignition, as shown in Fig. 6. High reactivity brought by the high molar concentration of the gas mixture produces robust thermal expansion and elevates the flame velocity above 300 m/s in a short time. At this speed, the role of the gas compression is no longer negligible and moderates the flame acceleration, as concluded by Bychkov el al. [28]. Later, pressure waves from nearby flame front enhance the suppression effect on the movement of the unburnt gas, this effect grows stronger as the interfacing flames getting closer,which finally forces the flame to decelerate.It is interesting that the left-spreading flame from Ignition 3 undergoes a short period of reverse propagation,which is mainly due to the stronger compression effect from the opposite flame.

Fig.7. Temperature(red lines)and pressure(black lines)distribution representing the flame front and the pressure peak formation of the right-spreading flame from Ignition 3. Time instants are from t1=0.1 ms to tend=0.42 ms, Δt=0.04 ms.

The profiles of the temperature and pressure distribution shown in Fig. 7 represent the flame acceleration process which spreads from ignition point 3 towards the base of the projectile. As the flame propagates forward, thermal expansion of burnt gas produces pressure waves continuously. Compression heating of the pressure waves creates a preheat zone in the unburnt region ahead of the flame, as shown in Fig. 6. The higher species concentration and temperature make the reaction and chemical heat release stronger, which helps the following pressure waves to rise more.The positive feedback of the shock-flame interaction makes a contribution to the development of a flame with higher speed.

Before detonation transition, a pocket of fully compressed and heated unreacted zone is formed ahead of the flame.As a result,the flame gets fed with a fresh mixture of higher density and reactivity,and hence generates a higher pressure wave and pushes the unburnt zone ahead of the flame to a higher velocity. Both the pressure peak and the flame velocity start to grow dramatically.Detonation transition happens from 0.4423 ms to 0.4435 ms aroundxDDT=659.6 mm,the pressure elevates from 196 MPa to an overdriven detonation value of 1063 MPa. The overdriven detonation decays as it propagates towards the projectile base,but due to the short propagation distance, the pressure is still above the C-J value when the detonation wave reaches the wall,as Fig.8 shows.

4.3. Comparison of the detonation ballistic performance against the normal ballistic cycle

It is noticed that detonation is triggered on the longest flame propagation path in the 3-point ignition case. Hence, normal ballistics which have no detonation phenomenon are simulated by rearranging shorter ignition spacing to study the influence of the detonation phenomenon on the ballistic performance.

Fig. 9 shows the flame velocity history of 5-point ignition and hydrogen-oxygen molar ratio of 6:1 until the flames are annihilated by collision with the flames propagating in the opposite direction or the end-wall.It can be seen that velocity evolution of the flame spreading from ignition 1 to the left end-wall stays similar to the 3-point ignition case.While the duration of the fast deflagration stage of the flames propagating in the opposite direction becomes shorter and the deceleration stage is shifted to an earlier point,as a result of the ignition spacing decreasing from 150 mm to 120 mm.The major difference occurs on the rightmost spreading flame,where the propagation distance is shortened from 350 mm to 170 mm. Flame propagation starts decelerating obviously from 0.42 ms when the flame front is about 12 mm from the bullet. As the flame propagates on this end-gas region, reflection waves become comparable to the new generated pressure waves on the flame front and hence inhibit the movement of the unburnt gas.During this period, although the average pressure of the unburnt gas rises quickly,it is difficult to form synchronous propagation and continuous enhancement between the flame and pressure wave as a result of the low flame velocity, hence, detonation is avoided.

Fig. 8. Attenuation profiles of the overdriven detonation pressure from t1=442 μs to tend=464 μs,Δt=2 μs, dash line is the C-J pressure calculated with the local gas state at t=442 μs.

Fig.9. Flame velocity evolution of 5-point ignition and MH2:MO2=6,the absent flame velocity curves of other propagating directions are similar to the shortest three curves.

Fig.10. Pressure records on the projectile base of the detonation occasion against the normal combustions.

Fig.11. Velocity records of the projectile of the detonation occasion against the normal combustions.

Comparisons of the projectile base pressure and projectile velocity are shown in Fig.10 and Fig.11.For the detonation occasion,the projectile moves immediately as the detonation wave reaches its base, the instantaneous pressure jumps to 721 MPa, violent impact causes an elevation on the initial acceleration. The detonation wave then reflects and attenuates in the chamber, forming low-frequency high-intensity pressure oscillations. Consequently,there can be seen fluctuations on the velocity curve.It takes 3.73 ms for the projectile to accelerate out of the barrel, and the muzzle velocity is 1859 m/s.

For the case of 5-point ignition and MH2:MO2=6,the projectile starts moving at 0.66 ms when the pressure on the projectile reaches 172 MPa.The unburnt gas then passes through the converging chamber throat and burns in the gun barrel,increasing the pressure on the projectile to the peak value of 277 MPa. The projectile accelerates more smoothly and it takes 4.67 ms for the projectile to move out of the barrel. Dilution of the hydrogen-oxygen mixture will reduce both the peak pressure and the muzzle velocity, i.e.changing hydrogen-oxygen ratio to 8:1 will reduce the peak pressure by 17%-238 MPa and decrease the muzzle velocity from 1738 m/s to 1622 m/s.

4.4. Summary and discussion

Generally, multi-point flame propagation in CLGG can be separated into three distinctive stages.The first stage is an exponential acceleration in the early stage after ignition,which flames from all ignition kernels will experience. During this stage, the flame velocity is mainly controlled by the velocity of the unburnt gas ahead of the flame surface and the density ratio of the unburnt to the burnt.Approximately,the propagating velocity can be written as an exponential function of the flow time, which isuflame=34.13 exp（48.694t/ms） for MH2:MO2=6. Alternative formation can be written as a detailed expression proposed by Bychkov et al. [28], i.e.where, R, B,σ are the laminar falme velocity of the initial unburnt gas mixture, chamber diameter, numerical factor and scaled acceleration rate.

The second stage is a fast deflagration stage, the transition occurs when a new balance is established between the thermal expansion and the resistance effect due to the gas compression at high flame speed. During this stage, the free expansion of the unburnt gas is also suppressed by the oncoming pressure waves from the opposite flame, hence no significant acceleration remains. As the flames get close at the end of this stage, pressure waves from the opposite flame become comparable to the new pressure waves generated by the flame front, consequently the flame propagation steps into the third fast deceleration stage. This phenomenon is more significant in the flame propagation process with delayed multipoint ignition [29], the strong compression waves from an earlier ignited flame will impede the expanding of the delayed one,forming a stagnation or reverse propagation flame.

Comparing the velocity curves of the leftmost flames in the 3-point and 5-point ignition cases, which have equal propagation distance and similar profiles, it can be concluded that adding ignition points on the other side of the ignition origin has a small effect on the flame acceleration. Flame velocity profiles of the middle ignition points indicate that multipoint flame propagation is highly concerned with the spacing between ignition points; a smaller spacing will lead to a shorter fast deflagration duration and hence the maximum flame velocity decreases. Since detonation transition requires the establishment of the continuous coupling enhancement between the shockwave and the reaction zone, a sufficiently long distance for the flame to accelerate to high velocity is necessary.From this point of view,adding more ignition points to achieve shorter propagation distances is an effective means to avoid detonation hazard.

5. Conclusion

Numerical simulations of premixed hydrogen-oxygen flame in a smooth tube and CLGG have been performed by solving the 1-D fully compressible Navier-Stokes equations coupled with a detailed chemical reaction model using the WAF scheme HLLC Riemann solver with grid refinement method.

Three distinct stages are found in the multipoint ignited flame propagation: an initial stage of exponential acceleration, a middle stage of constant acceleration, and a last stage of fast deceleration or cumulative reacceleration until detonation. Gas compression causes moderation of the initial exponential flame acceleration while the impedance effect of the pressure waves from face to face propagating flames on the unburnt mixture movement is the main reason for the fast deceleration.For a deflagration flame with high velocity, the reaction rate and heat release, the pressure peak generated due to the thermal expansion, and the unburnt gas velocity, temperature and density form a positive feedback, pushing the fast deflagration flame to a supersonic state and transforms into an overdriven detonation through the continuous couplingenhancement between the synchronously propagating reaction wave and the pressure peak.

The occurrence of detonation phenomenon will improve about 7%performance on the muzzle velocity against the normal ballistic cycle, but it also generates a drastic impact shock with more than two times higher pressure peak and following high pressure oscillations. Proper ignition arrangement and dilution are helpful to avoid the detonation phenomenon, as well as its safety hazard to the CLGG facilities.

Acknowledgments

This work has been supported by the Foundation of National Key Laboratory of Shock Wave and Detonation Physics (Grant No.6142A0302020517).