Numerical investigation on combustion in muzzle flows using an inert gas labeling method

International Journal of Heat and Mass Transfer 101 (2016) 91–103

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Numerical investigation on combustion in muzzle flows using an inert gas labeling method Qiongyao Qin, Xiaobing Zhang ⇑ School of Energy and Power Engineering, Nanjing University of Science and Technology, Nanjing 210094, PR China

a r t i c l e

i n f o

Article history: Received 11 February 2016 Received in revised form 22 April 2016 Accepted 4 May 2016

Keywords: Muzzle flows Precursor gas Combustion Tracing technique Oxygen supply

a b s t r a c t The influence of the precursor flow on combustion in muzzle flows is investigated. The fourth-order Runge–Kutta method is employed to solve the classical interior ballistics model, providing velocity for the projectile when it accelerates along the barrel. An inert gas labeling method is proposed. An additional species, helium, is chosen as the label to tracing the precursor gas which fills the barrel before the projectile starts. A high-resolution upwind scheme, AUSM+ (Advection Upstream Splitting Method), and detailed reaction kinetics model are employed to solve the multispecies Navier–Stokes equations with finite rate chemistry. The precursor flow generated by the precursor gas driven out of the barrel ahead of the projectile is simulated. The development of muzzle flow with chemical reaction is simulated. It is demonstrated from the results that the secondary temperature rise in the intermediate region behind the Mach disk is attributed to combustion in this area. It is found that the core of the precursor gas supplies oxygen for combustion at 150 ls after the projectile base leaves the muzzle. Furthermore, despite the disrupted precursor flow, the precursor gas is still united and gradually diffuses into the propellant gas. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction This investigation is motivated by the interest in the combustion phenomenon for a gun-launched supersonic projectile and the interest in how to investigate a certain part which has the same physical and chemical properties with the surroundings. After the projectile base leaves the muzzle, propellant gas of high temperature and high pressure is released into the ambient air, generating an unsteady flow which contains a lot of undesired phenomenon such as blast waves, acoustic waves, pressure waves, electromagnetic radiation, muzzle flash and smoke [1]. In addition, the expulsion of a column of air driven out of the barrel ahead of the projectile generates a precursor flow field, which makes the muzzle flow more complex due to the precursor gas blast shock [2]. Extensive investigations of muzzle flows were carried out in the past decades, both experimentally and numerically [3–8]. Schmidt measured the structure of the flowfields formed about the muzzle of a small caliber rifle during the firing using a time-resolved, spark shadow-graph technique. In his experiment, the coupling and uncoupling between the initial gas blast field and the propellant jet flows were observed [9]. Based on the former experiment, he ⇑ Corresponding author. E-mail addresses: [email protected] (Q. Qin), zhangxb680504@163. com (X. Zhang). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.05.009 0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.

investigated the flow at the muzzle of a gun during launch of finstabilized projectiles both experimentally and analytically. The effect of launch gas dynamic loadings upon projectile motion was calculated and was shown to compare favorably with experiment [10]. Taylor formulated a model to describe the details of the flowfield produced by the firing of a gun or mortar, and numerical results agreed reasonably well with the existing experimental data. Besides, the development of precursor flow was found out in the same manner with propellant gas flow [11]. Subsequently, a lot of research was conducted on the precursor flow and its influence on propellant gas flow. A temporal new computational technique was adopted by Moretti to solve the Navier–Stokes equations which describe the precursor flow. By comparing with experimental results, the results of their numerical tests show that the basic features of their technique, viz. grid definition, integration scheme and shock fitting, do indeed provide very high accuracy [12]. The Arbitrary Lagrangian–Eulerian (ALE) of Euler equations was solved by Jiang using the AUSMDV scheme (a mixture scheme of two forms of the AUSM), and the prominent characteristics including the propagation of first and second blast waves, the generation of bow shock wave and moving of the projectile, etc. were discussed in detail based on the numerical results [13]. Especially, interactions between the precursor flow field and propellant gas exhaust were observed experimentally by Schmidt [14]. The influence of

92

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

Nomenclature b c D e f k l L m M Mw p q R ^ R S t T u u

v v0 v 00

X Z

burning coefficient of propellant, m s
1 pa
1 mass fraction diffusion coefficient half of the web-thickness of propellant, m impetus, J kg
1 rate constant of reaction displacement, m length, m mass, kg symbol of species molecular weight pressure, Pa capacity of heat transmission universal gas constant, 8.314 J mol
1 K
1 Arrhenius molar rate of generation area, m2 time, s temperature, K velocity vector, m s
1 velocity, m s
1 velocity, m s
1 stoichiometric coefficient for reactant stoichiometric coefficient for product molar fraction relative thickness of the burnt

the precursor flow field was analyzed by comparing data taken at fully ambient conditions with that acquired in firings from an evacuated gun tube, and it was concluded that precursor flow must be taken into consideration when muzzle flow is studied. In terms of combustion in muzzle flow, Klingenberg conducted experiments of combustion phenomenon associated with the flow of hot propellant gases [15–17]. In these experiments, spatial and temporal distributions of temperature, velocity and pressure of a rifle of a caliber 7.62 mm were measured, and the muzzle flash is separated into three main luminous regions in space and time, i.e., the primary flash, intermediate flash and secondary flash. Schmidt examined the influence of both transient and three-dimensional flow on the probability of ignition in gun plumes [18]. Klingenberg investigated the influence of turbulent mixing on combustion of muzzle flows. Conclusions were drawn that oxygen from the entrained air is transported to the core of the flow in the vicinity of the Mach disk, i.e., the so-called ‘‘intermediate” flash, or ‘‘intermediate” flow region and combustion reactions take place there. A more vital conclusion was drawn that the ignition sequence leading to secondary muzzle flash is governed by an ignition temperature whose value is determined by two energy sources: the well-known shock heating and the combustion reactions in the intermediate flash [19]. Zhuo simulated the muzzle flows with base bleed projectile based on dynamic overlapping grids. In this simulation, a chemical reaction kinetic model was adopted to describe the chemical non-equilibrium flow [20]. As concluded by Schmidt, precursor flow must be taken into consideration when muzzle flow is studied [14]. Previous studies analyzed the influence of the precursor flow based on comparing data acquired from cases with and without a precursor flow field [14], or directly simulating the muzzle flow coupled with the precursor flow [12,13]. But, to date, there is no published research describing the detailed process of the precursor flow interacting with the propellant gas flow, especially the oxygen contribution of the precursor gas to the intermediate flash which directly affects

Subscript i j r x y

species species reactions in x-direction in y-direction

Superscript n burning exponent Greeks w

v n

g u x a s q k

l D

mass fraction of the burnt form characteristic quantity form characteristic quantity form characteristic quantity minor work coefficient charge weight, kg covolume, m3 kg
1 viscous stress, N m
2 density, kg m
3 thermal conductivity, W m
1 K
1 molecular viscosity coefficient, N s m
2 charge density, kg m
3

the ignition of secondary combustion [19]. From Fig. 1 [14], it is obvious that it is hard to locate the precursor gas, which makes it impossible to analyze the motion of the precursor gas and the entrainment of oxygen. To investigate the influence of the precursor flow on combustion in muzzle flow, especially, the oxygen contribution of the precursor gas, the entire process from the start of the projectile in the barrel to the complete development of the muzzle flow, shown in Fig. 2, is simulated. In order to distinguish the precursor gas from the surrounding environment, and to investigate the oxygen contribution of the precursor gas as well as the interaction between the precursor flow and the propellant gas flow, an inert gas labeling method inspired by the research of Grynkiewicz [21] is proposed. In his research, fluorescent indicators were used for biochemical studies of the physiological role of cytosolic free Ca2+. In this investigation, helium is chosen as the inert gas label because of its chemical inertness to trace the motion of precursor gas, which is similar to that of the fluorescent indicators tracing Ca2+. With the help of helium, it is easy to directly trace the motion of precursor gas, interaction between the precursor flow and the propellant gas flow, as well as the entrainment of oxygen. To take

Fig. 1. Experimental shadowgraph of muzzle flow with precursor flow.

93

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

Fig. 2. The entire process from the start of the projectile in the barrel to the complete development of the muzzle flow: (a) the start of the projectile; (b) the formation of the precursor flow; (c) the formation of the propellant gas flow.

the precursor flow into consideration, the classical interior ballistics equations are solved by means of a fourth-order Runge–Kutta method [22]. The AUSM+ scheme and detailed reaction kinetics model are adopted to solve the chemical non-equilibrium Navier–Stokes equations. Besides, the k
e model is chosen as the turbulent model. This numerical simulation is performed using FLUENT (V. 14.5). It costs a month to finish the compute using 16 cores in parallel. Based on the results, the development of precursor flow, the development of muzzle flow, the influence of the precursor flow on combustion in muzzle flow are analyzed. 2. Mathematical model 2.1. Classical interior ballistics model In this simulation, the velocity of the moving projectile, or the moving boundary, is given by the classical interior ballistics model [22]. It consists of five equations: Form function The 7-perforated propellant, shown in Fig. 3(a), is used in this investigation, which follows the geometric burning law. The burning of the propellant is divided into three stages: the burning surface is increasing and the function is the first one bellows; when the burnt thickness equals the half of the web-thickness, the propellant splits into small grains with triangular section, shown in

Fig. 3(b), and the burning surface is decreasing; the propellant is burnt over. The functions of three stages is as follows

8 > vZð1 þ nZ þ gZ 2 Þ ðZ < 1Þ > <
 w ¼ vs ZZ 1 þ ns ZZ ð1 6 Z < Z k Þ k k > > : 1 ðZ P Z k Þ

ð1Þ

where w is mass fraction of the burnt, Z is relative thickness of the burnt. v, n and g are form characteristic quantities. vs and ns are the form characteristic quantities at the stage when the propellant surface starts to decrease. Burning law

dZ ¼ dt

(

ðZ < Z k Þ ðZ P Z k Þ

b n p e

0

where

e

is

half

ð2Þ

of

the

web-thickness

of

the

propellant

(e ¼ 0:000355), b is the burning coefficient (b ¼ 1:86  10
8 ), n is the burning exponent (n ¼ 0:84), p is the mean pressure in the barrel. Momentum equation

Spd
ðF q þ F f Þ ¼ um

dv dt

ð3Þ

where S is the sectional area of the barrel, pd is the pressure at the projectile base, F q is the resistance in front of the projectile, F f is the frictional resistance, m is the mass of the projectile, v is the velocity of the projectile, u is the minor work coefficient. Motion equation

v¼

dl dt

ð4Þ

Energy equation

h Spðl þ lw Þ ¼ f xw
umv 2 2

ð5Þ

where f and x are the impetus and the charge weight respectively, h ¼ c
1 and c is the adiabatic exponent, and lw is the diameter shrunk length of the free volume and is given by

"

Fig. 3. The cross-sectional view of the 7-perforated propellant: (a) the crosssectional view before burning; (b) the cross-sectional view at the moment when it’s about to split.

lw ¼ l0 1


D

qp


D a


1

qp

! #

w ; l0 ¼

V0 x ; D¼ S V0

ð6Þ

94

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

where l0 is the diameter shrunk length of the chamber, qp is the solid propellant density, D is charge density, a is the covolume, V 0 is the volume of the chamber. Relative variables are introduced to make the equations dimensionless, which are

l vj p l¼ ; t¼ ; p¼ ; l0 fD l0

v v¼ vl

ð7Þ

, v  are dimensionless quantities, where l, t, p of the projectile and is computed by

vl

v l is the limit velocity

sffiffiffiffiffiffiffiffiffiffiffi 2f x ¼ hum

ð8Þ

E ¼ ½q1 u; . . . ; qN u; qu2 þ p; quv ; uðqE þ pÞ

T

ð13Þ

F ¼ ½q1 v ; . . . ; qN v ; quv ; qv 2 þ p; v ðqE þ pÞ

ð14Þ

 T @c1 @cN Ev ¼ qD1 ; . . . ; qDN ; sxx ; sxy ; usxx þ v sxy þ qx @x @x

ð15Þ

 T @c1 @cN ; . . . ; qDN ; sxy ; syy ; usxy þ v syy þ qy Fv ¼ qD1 @y @y

ð16Þ

T

H¼


The dimensionless equations are as follows,

8 8
qffiffiffiffi > h > > > v 1 þ 2nZ þ 3gZ 2 pn ðZ < 1Þ > > 2B > > dw <
> qffiffiffiffi > > ¼ vs 1 þ 2n Z h > pn ð1 6 Z < Z k Þ > s Zk dt > > Zk 2B > > > > : > > 0 ðZ P Z k Þ > > ( qffiffiffiffi > > > h n > dZ p ðZ < Z Þ > k 2B > ¼ > < dt 0 ðZ P Z k Þ > > dl > > ¼v > > > dt > > > > > dv ¼ h p
ðF q þ F f Þul0 h > > d > 2f x > > dt 2 ! # " > > > > dp 1 1 dw 1 þ h > > > : dt ¼ ðl þ l Þ 1 þ D a
qp p dt
l þ l p v w w where

lw ¼ 1


D

qp


D a


1

!

qp

w; B ¼

S2 e 2 f xumb

Hv ¼

_ N ; 0; 0; 0T _ 1; . . . ; x S ¼ ½x

ð9Þ

2ð1
nÞ

The flow inside of the barrel and outside of the muzzle is assumed to be two-dimensional axisymmetric. There is no reaction when the projectile accelerates in the barrel, so the Navier–Stokes equations without source terms are adopted to simulate the formation of the precursor flow. After the projectile base leaves the muzzle, the unsteady Navier–Stokes equations with chemical reactions are employed to simulate the formation of the muzzle flow. The equation in integral from is as follows: X

@Q dX þ @t

@X

Z

ðf  nÞdS ¼

X

ðH þ Hv þ SÞdX

ð11Þ

where, f ¼ ðE
Ev Þi þ ðF
Fv Þj, n is normal vector of the area element, Q is vector of the conservation variables, E and F are vectors of the convection fluxes, Ev and Fv are vectors of their corresponding viscous fluxes, H and Hv are vectors of the source term caused by axial symmetry, S is vector of the chemical reaction source terms and they are given by

Q ¼ ½q1 ; . . . ; qN ; qu; qv ; qET

ð12Þ

ð19Þ

where sxx , sxy , syy and shh are the viscous stresses, qx and qy are the capacities of heat transmission, and the expressions of them are as follows:

sxx ¼ 2l



@u 2 2 @u @ v v
lr  u ¼ l 2

@x 3 3 @x @y y

ð20Þ

syy ¼ 2l



@v 2 2 @ v @u v

lr  u ¼ l 2
3 @y 3 @y @x y

ð21Þ



@u @ v þ @y @x

ð10Þ

2.2. Governing equations

Z

 T 1 @c @c qD1 1 ; . . . ; qDN N ; sxy ; syy
shh ; usxy þ v syy þ qy y @y @y

sxy ¼ l ðf DÞ 2

ð17Þ

ð18Þ

These equations are solved utilizing a fourth-order Runge–Kutta method, providing velocity of the projectile when it accelerates in the barrel. Correspondingly, the pressure of the flow field is gained by solving the governing equations and by integrating the pressure over the head of projectile, the air resistance is acquired. In short, the results of the classical interior model provide velocity for the projectile, and the results of governing equations provide the value of F p for Eq. (3). This coupling between the interior ballistics model and the flow model is implemented by the User-Defined Function (UDF).

Z

T 1 q v ; . . . ; qN v ; quv ; qv 2 ; v ðqE þ pÞ y 1

v

ð22Þ

2 3

shh ¼ 2l
lr  u y

ð23Þ

qx ¼ k

N X @T @ci þ q Di hi @x @x i¼1

ð24Þ

qy ¼ k

N X @T @ci þ q Di hi @y @y i¼1

ð25Þ

where qi ði ¼ 1 . . . NÞ is the density of species i, u and v are the velocity component in x-direction and y-direction, respectively. p is the pressure of the fluid, E is the total energy of the fluid of unit mass, ci is the mass fraction of species i which is calculated by ci ¼ qi =q, T is the temperature, Di is the diffusion coefficient of species i which is expressed by

Di ¼ ð1
X i Þ

, X Xj j–i

Dij

ð26Þ

where, Dij is the binary mass diffusion coefficient of component i in component j, X i is the molar fraction of species i, l and k are the molecular viscosity coefficient and the thermal conductivity of the _ i is the net generation rate of species i gas mixture, respectively. x from chemical reactions, which is computed as the sum of Arrhenius reaction sources over the N r reactions,

x_ i ¼ Mw;i

Nr X ^ i;r R

ð27Þ

r¼1

^ i;r is the Arrhewhere Mw;i is the molecular weight of species i and R nius molar rate of generation of species i in reaction r.

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

In this simulation, the species and the mixture are assumed to be thermally perfect and the equations of state are

pi ¼ R

qi M w;i

T; p ¼ R

q Mw

ð28Þ

T

where Mw;i and M w are the molecular weights of species i and the mixture, respectively. The mass concentration of the mixture is the sum of the mass concentration of all species, the molar concentration of the mixture is the sum of the molar concentration of all species and the pressure of the mixture is the sum of the pressure P P P of all species, i.e., q ¼ Ni¼1 qi ; n ¼ Ni¼1 ni ; p ¼ Ni¼1 pi , so the relation between M w;i and M w is

Mw ¼

N X X i M w;i :

ð29Þ

i¼1

2.3. Chemistry kinetic model The chemical mechanism used in this study is taken from the work of Gibeling and Nietubicz [23], and basically consists of 8 species (CO, CO2, H2O, H2, O2, H, OH, O), 2 inert species (N2, He) and 12 reactions, which are given in Table 1, as well as the parameters for Arrhenius expression. The general form of reaction r is as follows: N X

v

i¼1

0 i;r M i

kf ;r

N X

kb;r

i¼1



95

The backward rate constant, kb;r , is computed using the following expression:

kb;r ¼

kf ;r Kr

ð33Þ

where K r is the equilibrium constant for reaction r and computed by N X


ðm00
m0 Þ  DSr DHr patm i¼1 i;r i;r K r ¼ exp
R RT RT

ð34Þ

where patm denotes atmospheric pressure (101,523 Pa). The term within the exponential function represents the change in Gibbs free energy, and its components are computed as follows: N
S DSr X m00i;r
m0i;r i ¼ R R i¼1

ð35Þ

N
h DH r X m00i;r
m0i;r i ¼ RT RT i¼1

ð36Þ

where Si and hi are the entropy and enthalpy of the species i evaluated at temperature T and atmospheric pressure. These values are given in the properties of the mixture material. 3. Numerical methods

v

00 i;r M i

ð30Þ

where m is the stoichiometric coefficient for reactant i in reaction r, m00i;r is the stoichiometric coefficient for product i in reaction r, Mi is 0 i;r

the chemical symbol of species i, kf ; r is the forward rate constant of reaction r and kb;r is the backward rate constant of reaction r correspondingly. The reactions involved in the muzzle flow are regarded as reversible reaction, thus the molar rate of generation of species i is given by

^ i;r ¼ Cðm00
m0 Þ kf ;r R i;r i;r

N Y

g0j;r

½cj;r 

j¼1


kb;r

N Y

!

½cj;r 

m00j;r

ð31Þ

j¼1

where C is the effect of third-body, which is neglected in current study, cj;r is the molar concentration of species j in reaction r, g0j;r is the rate exponent for reactant species j in reaction r. The forward rate constant of reaction r, kf ; r , is given by the Arrhenius expression

kf ;r ¼ Ar T br e
Er =RT

ð32Þ

where Ar is the pre-exponential factor, br is the temperature exponent, Er is the active energy for reaction r and R is the universal gas constant. Table 1 CO–H2–O2 combustion mechanism. Detailed reaction

Ar [consistent units] br [dimensionless] Er ½J=Kmol

H + O2 = OH + O H2 + O = OH + H O + H2O = OH + OH OH + H2 = H2O + H O + H + M = OH + M O + O + M = O2 + M H + H + M = H2 + M H2O + M = OH + H + M O2 + H2 = OH + OH CO + OH = CO2 + H CO + O + M = CO2 + M CO + O2 = CO2 + O

1.2e+17 1.5e+7 1.5e+10 1.0e+8 1.0e+16 1.0e+17 9.7e+16 1.6e+17 7.94e+14 4.4e+6 5.3e+13 2.5e+12


0.91 2.0 1.14 1.6 0
1.0
0.6 0 0 1.5 0 0

0.691e+8 0.316e+8 0.722e+8 0.138e+8 0 0 0 4.78e+8 1.87e+8
0.031e+8
0.19e+8 2.0e+8

3.1. Spatial discretization To obtain the complete process of the muzzle flow, which means the generation, development and attenuation, the moving boundaries have to move a distance equal to dozens or even hundreds of barrel diameters. A structured grid and dynamic layer method are applied to deal with the large displacement condition. For an arbitrary control volume, the spatial discretization is as follows:

Xc

@Q þ tðf  nÞdS ¼ Xc ðH þ Hv þ SÞ @t @ X

ð37Þ

where Xc is the volume of the control volume, @ X is the enclosed surface of the control volume, H, Hv and S are the values of H, Hv and S at centroid respectively. Due to the regular edges, the double integral can be replaced by the sum of ðf 1:j
f 2;j Þ  nj  Sj , that is

tðf  nÞdS ¼

@X

Nj X ðf 1;j
f 2;j Þ  nj  Sj

ð38Þ

j¼1

where f 1;j and f 2;j are the inviscid and viscous flux respectively, Sj , nj are the area and unit normal vector of j th surface of the control volume. Thus the spatial discretization can be expressed as follows: J  @Q 1 X ¼
f
f 2;j  nj  Sj þ H þ Hv þ S @t Xc j¼1 1;j

N

ð39Þ

3.2. AUSM+ scheme In the process of spatial discretization, the quantity of f that crosses the boundary Sj needs to be computed. In the present investigation, the inviscid flux is computed by the AUSM+ scheme [24,25]. This scheme splits the inviscid flux, E and F, into a convection term and a pressure term. Take F for example:

F ¼ Fc þ Fp

ð40Þ

where

Fc ¼ M  a  Q

ð41Þ

96

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

Q ¼ ðq; qu; qv ; qHÞT

ð42Þ

Fp ¼ ð0; p; p; 0ÞT

ð43Þ

where H ¼ e þ p=q is the enthalpy. From the expressions above, it can be seen that the convection term contains the Mach number M, local sound velocity a and flow parameter Q , and the pressure term only contains pressure p. The flux across the boundary can be expressed as:

Fc1=2 ¼ M 1=2  a1=2  Q 1=2

ð44Þ

Fp1=2 ¼ ð0; p1=2 ; p1=2 ; 0ÞT

ð45Þ

The flow parameter Q is taken from the control volumes on both sides of the boundary according to the flow direction at the boundary:

Q 1=2 ¼

QL

M1=2 P 0

QR

M1=2 < 0

ð46Þ

Considering the Mach number M is impacted by the characteristic wave, the splitting function is used to compute the Mach number at the boundary: 

M ¼

81 < 2 ðM  jMjÞ


2  1Þ  M 2
1

:  1 ðM 4

2

1 8

jMj P 1 jMj < 1

ð47Þ

The pressure splitting function is 

p ¼

81 < 2 ð1  signðMÞÞ : 1 ðM  1Þ2 ð2  M Þ 

3 M 16

4


2 M2
1

jMj P 1 jMj < 1

ð48Þ

In addition, the method of computing local sound velocity a1=2 has effect on the computational accuracy. To guarantee the accuracy of capturing shock wave, the local sound velocity of the boundary is computed with the sound velocities of both sides considered. The expression is as follows:

( a1=2 ¼

ða1=2 Þ2 =uL

juL j P a

minð e aL ; e a R Þ juL j < a

ð49Þ

3.3. Time-operator splitting method In the simulation of non-equilibrium chemically reacting flow, the characteristic time scale of the chemical reactions (Dt chem ) is much smaller than that of the flow (Dtflow ) and the characteristic time scales of all elementary reactions are different from each other, which is the so-called stiff problem. Splitting the characteristic time scale of the chemical reactions by using a time-operator splitting method can simplify the problem. The time step Dtflow is split into small parts which have the same magnitude as Dt chem and become the time step for the ordinary differential equations of the chemical reactions. The governing equation is divided into two parts:

Z X

Z X

@Q dX þ @t @Q dX ¼ @t

Z

Z

@X

ðf  nÞdS ¼

X

ðH þ Hv ÞdX

ð50Þ

Z X

ðSÞdX

ð51Þ

Firstly, the governing equation is computed with the chemical reactions frozen, from which the flowfield parameters are gained. In this process, the time step is Dt flow . Secondly, the ordinary

differential equations of chemical reactions are computed with the time step Dt chem . This integration should not stop until the total integration time equals Dt flow . The linear stability theory used in the FLUENT shows that the density-based implicit formulation adopted in this simulation is unconditionally stable, so the CFL (Courant number) is chosen 5.0. For supersonic flow, the time step Dt for every cell is computed from the CFL (Courant–Friedrichs–Lewy) condition

Dt ¼ CFL

X

ðjuj þ aÞSx þ ðjv j þ aÞSy

ð52Þ

where X is the area of the grid cell, u and v are the velocity components of the fluid, a is the sound velocity, Sx and Sy are the projection area of grid cell in x and y direction, respectively. The entire time step for solving the flow parts is the minimum time step for all grid cells. 4. Computational model The computational domain and the corresponding boundary conditions are shown in Fig. 4. The length of the barrel is L2 ¼ 1:7 m (the travel of the projectile is 1.502 m and the diameter shrunk length of the chamber is 0.198 m), the radius of the barrel is R2 ¼ 0:015 m (the radius of the projectile is 0:015 m as well); the length of the domain outside of the muzzle is L1 ¼ 1:2 m, the radius of this domain is R1 ¼ 0:8 m. The breech (indicated by ‘‘1”), the projectile base, the projectile head and the muzzle (indicated by ‘‘5”) are ‘‘wall”. The boundary of the outside domain is ‘‘pressure outlet”, including the boundary indicated by ‘‘4”. Attention needs to be paid to the boundary ‘‘interface”. Before specifying the ‘‘interface” boundary, it’s necessary to introduce the dynamic mesh used in this simulation. The dynamic layering method is employed to dealing with the moving boundary. In this method, the projectile base and head are specified as rigid body, of which the velocity is computed from interior ballistics when the projectile accelerates in the barrel and from the Newton’s Second Law with the pressure around the base and head after the base exits the muzzle. Besides, the mesh of the regions, ‘‘2” and ‘‘3”, is specified as rigid body and moves at the same speed with the projectile, which means the mesh splits at the breech and merges at boundary ‘‘4”. It can be seen from the introduction above that the regions on both sides of the dash need data exchange which is implemented by the ‘‘interface”. As mentioned above, the mesh used in this simulation is structured mesh which is shown in Fig. 5. Besides, unstructured mesh is adopted to deal with the curved surface of projectile head, which is shown in the magnified area in Fig. 5. The size of the structured mesh cell is 2 mm  2 mm, the size of the unstructured mesh cell ranges from 1 mm to 2 mm as shown in Fig. 5, and the total cells are 213,500. The size of the mesh comes from a grid sensitivity investigation in which a series of mesh size is investigated and the mesh introduced above shows great mesh independence. The simulation consists of two stages: the first stage is the acceleration of the projectile in the barrel and at this stage the precursor flow is formed. Before the start of the projectile, the pressure of the area ahead of the projectile is 101,325 Pa, the temperature is 300 K, the velocity is zero, the mass fractions of the species of the area ahead of the projectile in the barrel are O2 (0.233), N2 (0.667), He (0.1) and the mass fractions of the species of the area outside the muzzle are O2 (0.233), N2 (0.767). The velocity of the projectile is computed from interior ballistics model, as mentioned above, so the initial conditions of the area behind the projectile are set appropriate for computing. When the projectile is about to exit the muzzle, its velocity is up to 894 m/s and it comes to the next stage, the development of the muzzle flow.

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

97

Fig. 4. Schematic of the computational domain.

Fig. 5. Mesh setup.

Table 2 Mass fraction of species of propellant gas.

5. Inert gas labeling method

Species

CO

H2

CO2

H2O

N2

Propellant gas

0.5138

0.0157

0.2153

0.12993

0.1259

Before the exit of the base, the pressure of the propellant gas behind the base is set 50 MPa, the temperature is 1800 K and the velocity is 894 m/s. The components of the propellant gas are given in Table 2 [20].

In this investigation, a new method, inert gas labeling method, is proposed to trace the precursor gas. By this method, the motion of the precursor gas, the entrainment of oxygen and the interaction between the precursor gas and the propellant gas are traced by the label. The mass fraction of the helium in this investigation is 0.1. Due to that the helium is added into precursor gas, it is necessary to test whether the helium will influence the flow and the reactions. Therefore, a control simulation which does not contain

Fig. 6. Mach number, pressure and temperature contours of the labeled and unlabeled flow field at t ¼ 200 ls after the base leaves the muzzle: the upper half is the labeled and the lower half is the unlabeled.

98

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

Fig. 7. Pressure contours at difference time from 2.6 ms to 3.0 ms.

Fig. 8. Mach number contours at difference time from 2.6 ms to 3.0 ms.

Fig. 9. Mass fraction of helium contours at difference time from 2.6 ms to 3.0 ms.

helium is designed to investigate the influence. In Fig. 6(a)–(c), the Mach number, pressure and temperature contours of the labeled and unlabeled flow field are compared. From these figures, no obvious difference is shown in the Mach number, pressure and temperature contours, so conclusion can be drawn that the helium does not affect the flow and the reactions.

6. Results and discussion 6.1. The precursor flow field The projectile, of which the velocity is obtained by solving the classical interior ballistics model, accelerates inside of the barrel

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

99

Fig. 10. Mach number contours of muzzle flow.

Fig. 11. CO2 mass fraction contours of muzzle flow.

Fig. 12. Temperature contours of muzzle flow.

and continuously compresses the air column ahead of the projectile, generating a planar shock. A precursor flow field generates when the compressed air column behind the plane shock exits the muzzle. Fig. 7(a)–(c) shows the pressure contours at different time and the pressure shown in Fig. 7(c) at 3.0 ms after the start of the projectile is at the same magnitude with the pressure of ambient air. Fig. 8(a)–(c) shows the Mach number contours at different time. Especially, a structure similar to the typical feature of muzzle flow, i.e., Mach disk, is observed in Fig. 8(b), which is

subsequently destroyed by the projectile in Fig. 8(c). In Fig. 9(a)–(c), the precursor gas labeled by helium is traced, which indicates the state of oxygen in precursor gas. 6.2. Development of muzzle flow The exit of the high-pressure and high-temperature propellant gas out of the muzzle generates a highly under-expanded supersonic jet which interacts with the ambient air and forms the

100

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

Fig. 13. Comparison between mass fraction contours at t ¼ 400 ls.

Fig. 14. Temperature and mass fraction distribution along y ¼ 30 mm at t ¼ 400 ls.

Fig. 15. Comparison between He and CO2, O2, temperature.

muzzle flow. Fig. 10(a)–(c) shows the development of muzzle flow from 200 ls to 900 ls. Fig. 10(a) shows the Mach number contour of muzzle flow at 200 ls after the exit of projectile base. The flow is restrained by the projectile, forcing the highly under-expanded supersonic jet to form a bow shock which is indicated by ‘A’. As the propellant gas continuously jets into the muzzle flow and the

projectile is moving forward, the flow gradually frees from the restraint and eventually generates an unrestrained triple point and a fully developed inner shock disk (Mach disk) which are indicated in Fig. 10(b) by ‘B’ and ‘C’, respectively. Fig. 10(c) shows the subsequent development of muzzle flow. Fig. 11(a)–(c) gives the CO2 mass fraction contours which can show the status of propellant gas.

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

Fig. 16. Temperature and mass fraction distribution along y ¼ 30 mm at t ¼ 400 ls.

6.3. Combustion in muzzle flow Fig. 12(a)–(c) shows the temperature contours at t ¼ 200 ls, t ¼ 400 ls and t ¼ 900 ls, respectively. In these contours, the yellow and red region represent high temperature, which are indicated by ‘D’, ‘E’ and ‘F’, respectively. It can be seen clearly in Fig. 12(b) that the temperature in the flow decreases just downstream of the muzzle and increases when the flow crosses the

101

Mach disk. Subsequently, the flow undergoes a secondary increase in temperature downstream of the Mach disk. Although in previous investigations [15–17], the reactions taking place in the muzzle were confirmed, here a brief comparison is necessary to clarify that the secondary temperature increase is induced by chemical reactions. To make the analysis clear, the reaction region surrounded by a black pane in Fig. 12(b) is magnified. The combustion region is the core of intermediate flame, as mentioned by Klingenberg [19]. The magnified area is shown in Fig. 13(a)–(c). In Fig. 13(a)–(c), the comparisons between temperature and CO2 mass fraction, CO and CO2 mass fraction as well as CO2 and O2 mass fraction in the magnified area are made to demonstrate chemical reactions induce the secondary temperature increase. In Fig. 13(a), the upper contours are temperature and the lower CO2 mass fraction. Three comparative points (indicated by ‘‘1”, ‘‘2” and ‘‘3”) are chosen to compare the locations of the temperature increase and CO2 mass fraction increase. The black dash lines are auxiliary lines that help to confirm that the upper points and the lower ones are at the same position in flow field. From Fig. 13(a), it can be seen that where the temperature increases, the CO2 mass fraction increases. In Fig. 13(b), the mass fraction of CO2 shown in the lower half is compared with that of CO in the upper half utilizing the same method adopted above. The three upper points locate the surface of the CO region while the others locate the CO2 increase region, which indicates that CO2 increases at the surface of the CO region. In Fig. 13(c), the mass fraction of

Fig. 17. Mass fraction distribution along y ¼ 30 mm at different time.

102

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103

(a) t = 200 µs

(b) t = 400 µs

(c) t = 900 µs

Fig. 18. Precursor gas labeled by He.

CO, the upper one, and the mass fraction of O2, the lower one, are compared. Obviously, the surface of the CO region and that of the O2 region are in contact with each other. To confirm that the secondary temperature increase is induced by the reactions, temperature and mass fraction distributions along y ¼ 30 mm from x ¼ 1:96 m to x ¼ 2:1 m (indicated by the black solid line in Fig. 13(b)) at t ¼ 400 ls are shown in Fig. 14 (a) and (b). The mass fractions of OH and O, shown in Fig. 14(a), reflect the local reaction rate. O and OH appear between x ¼ 2:01 m and x ¼ 2:03 m, which indicates that reactions occur in this period along y ¼ 30 mm. Obviously, this is where the secondary temperature increase happens. That the temperature increase happens before x ¼ 2:01 m is due to heat conduction. Fig. 14(b) shows temperature and the mass fractions of CO2, CO, O2. The decreasing to zero between x ¼ 2:01 m and x ¼ 2:03 m of CO and O2 indicates the surface of the CO and O2 region. The cross of the two lines indicates mixing of CO and O2.

6.4. Contribution of precursor flow Fig. 15(a)–(c) shows comparisons between CO2 and He, temperature and He, O2 and He. As utilized above, the comparative method is adopted to compare three chosen points (‘‘1”, ‘‘2” and ‘‘3”). By comparison, the region where the CO2 and temperature increase and where the surface of O2 is located is the region filled with He, which indicates that oxygen for combustion is supplied through the region filled with He. Besides, Fig. 16 shows the mass fraction of He and O2 distributions along y ¼ 30 mm at t ¼ 400 ls. The distribution of O2 and He also demonstrates that oxygen is supplied through the He region. Fig. 17(a)–(d) shows a series of mass fraction of He and O2 along y ¼ 30 mm at t ¼ 100 ls, t ¼ 150 ls, t ¼ 200 ls and t ¼ 400 ls after the projectile base leaves the muzzle. In these figures, the dash and dot lines are auxiliary lines which help to confirm the relative location of the helium peak and oxygen diffusion front. The helium peak is the maximum of the mass fraction of He. The diffusion front is the location where diffusion of O2 occurs. With the help of the auxiliary lines, it can be seen in Fig. 17(a) that O2 diffusion occurs at x ¼ 1:850 m while the He peak locates at x ¼ 1:8625 m. The peak represents the core of the precursor gas, so a conclusion can be drawn that the core of the precursor gas does not supply oxygen for combustion at t ¼ 100 ls. In Fig. 17 (b), O2 diffuses at x ¼ 1:8975 m, and that is where the He peak is located. This indicates that oxygen in the core of precursor gas begins to diffuse into the combustion area at t ¼ 150 ls. In Fig. 17(c) and (d), it is obvious that diffusion occurs behind the

peak, which indicates that the gas behind the core of the precursor supplies oxygen through the helium. The comparisons between the CO mass fraction and the He mass fraction at different time are shown in Fig. 18(a)–(c). Schmidt has concluded that the propellant gas flowfield is highly energetic and rapidly expands through the precursor flowfield, effectively destroying it [9]. Although the typical structure, Mach disk, is destroyed by the projectile, it can be seen from the figures that the precursor gas is in contact with the propellant gas and restrains the propellant gas. 7. Conclusions The process from the start of the projectile in the barrel to the complete development of the muzzle flow is simulated. The development of precursor flow and muzzle flow is obtained. The motion of precursor gas is traced using an inert gas labeling method. The contribution of the precursor flow on combustion in the muzzle flow is analyzed. By comparing the location of CO2, CO, O2, He and temperature rise, the contribution of combustion on the temperature increase is demonstrated and the oxygen contribution of the precursor gas on combustion is also demonstrated. The change of the relative location of O2 diffusion front and the He peak is observed. Based on the discussion above, a conclusion is drawn that, the core of the precursor gas contributes oxygen for combustion at t ¼ 150 ls. The motion of the precursor gas is obtained. The precursor gas is in contact with the propellant gas behind the Mach disk and restrains the propellant gas. Acknowledgement This work is supported by the Natural Science Foundation of Jiangsu Province (Grant No. BK20131348) and Key Laboratory Foundation of the People’s Republic of China (Grant No. 9140C300206120C30110), the National Natural Science Foundation of China (Grant No. 11502114), China Postdoctoral Science Foundation funded project (Grant No. 2015M581797). References [1] G. Klingenerg, Gun muzzle blast and flash, Propellants Explos. Pyrotech. 14 (1989) 57–68. [2] J.I. Erdos, P.D. Del Guidice, Calculation of muzzle blast flowfields, AIAA J. 13 (8) (1975) 1048–1055. [3] R.A. Carson, O. Sahni, Numerical investigation of propellant leak methods in large-caliber cannons for blast overpressure attenuation, Shock Waves 24 (2014) 625–638.

Q. Qin, X. Zhang / International Journal of Heat and Mass Transfer 101 (2016) 91–103 [4] M.C. Taylor, T.L. Laber, B.P. Epstein, D.S. Zamzow, D.P. Baldwin, The effect of firearm muzzle gases on the backspatter of blood, Int. J. Legal Med. 125 (2011) 617–628. [5] X.H. Jiang, B.C. Fan, H.Z. Li, Numerical investigations on dynamic process of muzzle flow, Appl. Math. Mech. 29 (3) (2008) 351–360. [6] L.A. Florio, Effect of vent opening area and arrangement on gas flow field as gas propelled cylinder exits a flow tube, Meccanica 45 (2010) 475–501. [7] T. Mizukaki, Visualization and force measurement of high-temperature, supersonic impulse jet impinging on baffle plate, J. Vis. 10 (2) (2007) 227–235. [8] H. Rehman, S.H. Hwang, B. Fajar, H. Chung, H. Jeong, Analysis and attenuation of impulsive sound pressure in large caliber weapon during muzzle blast, J. Mech. Sci. Technol. 25 (10) (2011) 2601–2606. [9] E.M. Schmidt, D.D. Shear, Optical measurements of muzzle blast, AIAA J. 13 (8) (1975) 1086–1091. [10] E.M. Schmidt, K.S. Fansler, D.D. Shear, Trajectory perturbations of fin-stabilized projectiles due to muzzle blast, J. Spacecraft 14 (6) (1977) 339–344. [11] T.D. Taylor, T.C. Lin, Numerical model for muzzle blast flowfields, AIAA J. 19 (3) (1981) 346–349. [12] G. Moretti, A numerical analysis of muzzle blast precursor flow, Comput. Fluids 10 (1) (1982) 51–86. [13] X.H. Jiang, Z.H. Chen, B.C. Fan, H.Z. Li, Numerical simulation of blast flow fields induced by a high-speed projectile, Shock Waves 18 (2008) 205–212. [14] E.M. Schmidt, R.E. Gordnier, K.S. Fansler, Interaction of gun exhaust flowfields, AIAA J. 22 (4) (1984) 516–517. [15] G. Klingenberg, H. Mach, Investigation of combustion phenomenon associated with the flow of hot propellant gases-I: spectroscopic temperature measurements inside the muzzle flow of a rifle, Combust. Flame 27 (1976) 163–176.

103

[16] G. Klingenberg, G.A. Schröder, Investigation of combustion phenomenon associated with the flow of hot propellant gases-II: gas velocity measurements by laser-induced gas breakdown, Combust. Flame 27 (1976) 177–187. [17] G. Klingenberg, Investigation of combustion phenomenon associated with the flow of hot propellant gases-III: experimental survey of the formation and decay of muzzle flow fields and of pressure measurements, Combust. Flame 29 (1977) 289–309. [18] E.M. Schmidt, Ignition in gun exhaust plumes, Tech. Notes (1984) 441–443. [19] G. Klingenberg, J.M. Heimerl, Combustion following turbulent mixing in muzzle flows, Combust. Flame 68 (1987) 167–175. [20] C.F. Zhuo, X.S. Wu, F. Feng, Q. Liu, H. Ma, Numerical simulation of the muzzle flows with base bleed projectile based on dynamic overlapped grids, Comput. Fluids 105 (2014) 307–320. [21] G. Grynkiewicz, M. Poenie, R.Y. Tsien, A new generation of Ca2+ indicators with greatly improved fluorescence properties, J. Biol. Chem. 260 (6) (1985) 3440– 3450. [22] X.B. Zhang, Z.M. Jin, Interior Ballistics of Guns, Beijing Institute of Technology Press, Beijing, 2014. [23] H.J. Gibeling, C.J. Nietubicz, Navier–Stokes computations for a reacting, M864 base bleed projectile, AIAA Paper 93-0504, in: AIAA 31st Aerospace Sciences Meeting, Reno, NV, 1993. [24] M.S. Liou, A Sequel to AUSM: AUSM+, J. Comput. Phys. 129 (2) (1996) 364–382. [25] C. Shen, F.X. Sun, X.L. Xia, Analysis on transient conjugate heat transfer in gapcavity-gap structure heated by high speed airflow, Int. J. Heat Mass Transfer 67 (2013) 1030–1038.