- Email: [email protected]
A numerical study of a detonation wave in detonation transmission tubing
Mathematical and Computer Modelling 44 (2006) 717–734 www.elsevier.com/locate/mcm
A numerical study of a detonation wave in detonation transmission tubing D.K.L. Tsang ∗ School of Mechanical, Aerospace and Civil Engineering, The University of Manchester, PO Box 88, Sackville Street, Manchester M60 1QD, UK Received 12 July 2004; received in revised form 15 December 2005; accepted 9 February 2006
Abstract In this paper we formulate a one-dimensional detonation wave model for detonation transmission tubing. The reaction is described by a temperature-independent single-step kinetic equation and the resulting equations are a system of hyperbolic conservation laws. The initial boundary value problem is then solved by using the TVD MacCormack method. The effects of heat loss to the tube wall and of the tube expansions have also been considered. The result obtained shows that the heat loss and the tube expansion have a major effect on the pressure prediction. c 2006 Elsevier Ltd. All rights reserved.
Keywords: Detonation; Rarefaction wave; Detonation transmission tubing; Deflagration; HMX
1. Introduction Detonation transmission tubing (DTT) is a modern method for transmitting the detonation signal from a start device to the detonator delay element and hence to the main explosive. It is much safer than electric wiring, since conventional wire can be accidentally initiated by an electric discharge. The DTT consists of a tube formed from SURLYN, a copolymer of ethylene and methacrylic acid. The internal and external tube diameters are 1.3 mm and 3.0 mm, respectively. The interior surface is coated with a mixture of aluminium (Al) flakes and HMX. Core loadings are typically 16 mg/m (milligrams per meter length of the tube). About 92% of the mixture is the HMX and this is the main provider of energy to support the detonation wave. The rest is aluminium, which is mainly present to expedite initiation. Each HMX particle is approximately in the form of a cube with sides of length 40 micron. Each Al flake is in the form of a thin plate approximately 10 micron square and 0.1 micron thick. The mixture occupies about 1% of the internal tube volume. The DTT is normally initiated by a shock wave at one end. Associated with the shock wave are pressure, temperature and density changes. It is posited that the passage of the shock wave raises the temperature of the gas in the tube, and this subsequently causes the HMX to react. The temperature rise also heats the Al particles sufficiently for the oxide layer on the surface to soften and for rapid reaction between the exposed metal and oxygen to commence. There is interest at the present moment in trying to predict the effects of any accidental damage that may occur to the ∗ Tel.: +44 1612754506; fax: +44 1612754901.
E-mail address: [email protected]. c 2006 Elsevier Ltd. All rights reserved. 0895-7177/$ - see front matter doi:10.1016/j.mcm.2006.02.008
718
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
Nomenclature W ρal ρHMX ρp Dal DHMX Dp Kp C1 C2 C3 Cp 1H Rgas R L L1 T∗ T0 Ti p0 a0 b0 Ep νp
average molecular weight 22.99 g/mol density of aluminium 2700 kg/m3 density of HMX 1900 kg/m3 density of tube wall 920 kg/m3 diffusivity of aluminium 9.76 × 10−5 m2 /s diffusivity of HMX 1.00 × 10−7 m2 /s diffusivity of tube wall 1.03 × 10−7 m2 /s conductivity of tube wall 0.198 W/m K specific heat of HMX 1500 J/kg K specific heat of Al 903 J/kg K specific heat of gas 1000 J/kg K specific heat of tube wall 2100 J/kg K heat of reaction 5.66 × 106 J/kg gas constant 8.31 J/mol K perfect gas constant 361.65 J/kg K length of tube 2 m initial length of temperature variation 0.1 m cut-off temperature of HMX 450 K reference temperature 293 K ignition temperature 500 K reference pressure 1.01 × 105 Pa initial inner radius 6.5 × 10−4 m initial outer radius 1.5 × 10−3 m Young’s modulus of tube 1.80 × 108 Pa Poisson’s ratio of tube 0.49
tube during mining operations. For this purpose it would be useful to have a model of the operation of the tube and to understand the mechanism of the detonation process. Detonation wave is a combustion process initiated by a shock. Detonations were first investigated scientifically by Berthelot and Vielle [1,2], and by Mallard and Le Chatelier [3]. Hypotheses leading to mathematical predictions of the speed of propagation of the waves were given by Chapman [4] and Jouguet [5,6]. A good approximation of the usual detonation structure is a shock wave followed by a deflagration. This shock–deflagration structure of detonation was first proposed by Zeldovich [7], von Neumann [8] and Doering [9] and is referred to as the ZND wave structure. On the basic of numerical calculations of steady, planar detonation structure (for example, Oppenheim and Rosciszewski [10]) and of good experimental measurements performed mainly in the 1960s (for example, Edwards et al. [11]), it was widely believed that the ZND structure was representative of most real detonations. However, more recent studies established that this structure is usually unstable. In fact, these detonations generally possess a complex three-dimensional substructure involving transverse waves that produce a steady, planar wave only on average (see Lee [12]). In Abouseif and Toong [13], the mechanism of instability is identified. In their analysis, they demonstrated that the interaction between the irreversible temperature fluctuations and the finite reaction zone induces an oscillating energy source field, which then leads to shock perturbations and thereby detonation pressure oscillations. Detonation propagation in premixed substances is one of the basic aspects of combustion. The complexity of detonations is greatly increased if the pre-mixed system contains finely pulverized organic or metallic solids suspended in air or oxygen. Self-sustained detonations can be obtained in heterogeneous dust fuel-gas oxidizer mixtures, and the propagation mechanism is dominated by a transverse wave structure similar to those in homogeneous gas mixtures. These facts have been firmly established by the observations of Zhang and Gronig [14]. In their reactive corn starch particles/oxidizing gas system, they find that the momentum and temperature relaxations are dominant immediately
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
719
behind the shock front. There, the gas velocity and temperature jump at once to the new state corresponding to the shock strength, while the particle velocity and temperature increase continuously from the initial state. This causes the induction zone length to increase. Heterogeneous effects like pyrolysis, surface oxidation, etc, also cause the induction zone length to increase. Experimental observations for the DTT, provided to us by ICI Explosives Division, proffer that the detonation wave passes down the DTT with a typical speed of 2 km/s and a peak pressure of about 10–12 MPa. In the present work, we begin with the simplest hypothesis and formulate a one-dimensional wave model for the passage of the detonation wave along an unbroken length of DTT. The reactants are considered to be homogeneously distributed within the tube with the reaction rate described by a temperature-independent single-step kinetic equation. This results in a set of hyperbolic conservation equations that we solve numerically. The effect of heat loss to the tube wall and the effect of tube expansion have also been considered. The model showed good agreement with experimental observations. 2. Assumptions In our wave model, the tube wall interaction terms have been considered in the derivation. We now list and discuss the key assumptions made in the gas/particulars motion and in the tube wall of our one-dimensional model. The various time scales used in the discussion are derived in the Appendix. 2.1. Assumptions for gas/particulars motion The assumptions for the gas/particulars motion are as follows: 1. The gas flow is one-dimensional, so that all variables depend only on x, the distance down the transmission tube and on t, the time. In practice, when a detonation wave propagates in a tube, the front structure is intrinsically two- or threedimensional and non-steady. Also, the wall cooling effect, the wall expansion effect and the heterogeneous nature of the HMX prevents the problem from being one-dimensional. However the one-dimensional detonation wave is the simplest representative system, is significantly easier to solve, gives insight into the basic physics of the problem, and needs to be understood before investigations of any non-planar phenomena can properly begin. Therefore, we begin by considering the one-dimensional calculation. 2. The speed of the aluminium flake is the same as the gas speed, but the speed of the HMX is negligible, i.e. the HMX remains stationary at its original position. A typical relaxation time for the speed of the aluminium flake is 0.193 µs. We can see that this time scale is short compared to a reaction time scale of 40 µs. We can thus assume that the aluminium flake is instantaneously accelerated by the gas flow, and so always has the same speed as the gas itself. However, a typical speed relaxation time for the HMX particles is 27.1 µs. Clearly, the assumption that the HMX is stationary is not perfect. The particles would be accelerated to gas speed over a time comparable to the time of the chemical reaction. However, during the time that the conversion of HMX to gas is a maximum, the particle motion is small, so we would expect the HMX motion effect to be small. 3. The temperature of the HMX remains constant during the reaction, but the temperature of the aluminium flake is assumed to be the same as the gas temperature. The time scale to heat up an HMX particle by diffusion is 4.0 × 10−3 s. This is slow compared to the time scale for the reaction. We can conclude from this that the bulk of the particle does not have time to heat up by diffusion. Hence it does not react as a uniform whole, and must therefore react by conversion of the surface of HMX to gas. Consequently, the temperature of the unburnt HMX remains constant until reaction occurs. Conversely, the time scale taken to heat up an aluminium flake is 2.56×10−11 s. Now this value is very short compared with our reaction time scale, so we assume that the temperature of the aluminium flake is the same as that of the gas. 4. The rate of reaction Γ (defined as the rate of change of mass of HMX per unit volume of tube) is zero below a certain cut-off gas temperature T ∗ . Above that gas temperature, Γ is proportional to the surface area of HMX remaining. The kinetics for a slab of uniformly heated HMX has been investigated by [15], but these cannot be applicable here. We have already given reasons why the reaction is expected to be a surface reaction, rather than a volume
720
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
reaction, and for this reason we assume that the rate of reaction Γ is proportional to the surface area of HMX remaining. It is known that the HMX decomposes at 450 K, and we choose this as our cut-off temperature T ∗ . 5. Only a single chemical reaction takes place. The heat of this reaction is the combined heat of reaction of HMX decomposition and aluminium oxidization. The aluminium is treated as an inert phase which accumulates mass as the reaction proceeds. This assumption was made to reduce the number of chemical reactions that the eventual computer program would have to keep track of, and hence to produce a more tractable model. A number of simplifications follow from the assumption. The final state of chemical equilibrium can be established from thermodynamic considerations, and in the model it is assumed that the reaction progresses evenly towards this state in proportion to the amount of HMX that has decomposed. In the real process, oxygen is likely to combine preferentially with the aluminium present. Only when this has been completely converted to Al2 O3 would there be significant carbon monoxide and water products. In the same way, the oxygen initially present in the tube would probably combine very early on, rather that continuously over the lifetime of the HMX decomposition. Thus, the real chemical reaction would be multi-stage. 6. During reaction, 1 unit mass of HMX produces θ unit masses of oxygen and 1 − θ unit masses of inert gas, where θ is a constant. The θ unit masses of oxygen then react with the aluminium instantaneously, such that no free oxygen is present in the model. The chemical kinetics of HMX/Al system are very complex and unknown to us. In order to produce a tractable model, we assume only a temperature-independent single-step reaction and that the heat of this reaction is the combined heat of HMX decomposition and Al oxidization. We write the overall chemical formula per unit mole of HMX as follows: HMX + 0.954Al + 0.704N2 + 0.155O2 → 0.477Al2 O3 + 4CO + 2.879H2 O + 1.121H2 + 4.704N2 .
(1)
We treat the gases as a single gas product, taking the W molecular weight of the hypothetical gas to be the mole average of the final gas products, equal to 22.99 g/mol. The mass ratio θ is equal to 0.061. We can see that only a very small fraction of gas is absorbed by the aluminium flake. Most of the HMX becomes gas. 2.2. Assumptions for tube wall The assumptions for the tube wall are as follows: 1. The tube deformation is small. This assumption implies that the displacements are sufficiently small so that no distinction is needed between the coordinates of a particle before and after deformation, and that the displacement gradients are sufficiently small so that their products may be neglected. 2. The material behaves elastically at all times. This assumption implies that neither the temperature changes nor the stresses are too large. However, in our present study there is a huge difference between the internal pressure and external pressure. From our numerical results, we find that there is about a 14% increase in the inner tube radius at the peak pressure. So the stress changes are large and the deformation is not small. Nevertheless, in order to produce a tractable model, we have made this assumption. 3. The tube is displaced only radially at any locations. This is a plane strain assumption, and applies when the length of the body is large compared with the maximum cross-sectional dimension. In this study, we have a long tube with very small cross-sectional area, so the concept of plane strain may be used. 4. The thermal expansion is negligible. If we consider the plane strain problem, then one of the stress–strain relations and one of the strain–displacement relations are 1 e e εθθ = σθθ − νe σrr + αT (2) Ee ue εθθ = (3) r
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
721
where ε is the strain, σ e is the stress, u e is the radial deformation displacement, r is the radial coordinate, E p = 1.8 × 108 Pa is the Young’s modulus, ν p = 0.49 is the Poisson ratio, and α = 2.5 × 10−4 /K is the thermal expansion coefficient. From an order of magnitude analysis, we have [ p] + α [T ] = 0.07 + 0.75 [εθθ ] ∼ Ep where [ p] = 12 MPa is the peak pressure and [T ] = 3000 K is the peak temperature. We can see that the strain is dominated by the thermal expansion and that, at r = a0 , the deformation displacement [u] = 0.82a0 by using Eq. (3). However, such a large displacement is not observed. Presumably this is because the heated layer is above its melting point and undergoes molecular rearrangement, rather than elastic displacement. If this is correct, we would expect to see no thermal expansion at all but only physical expansion due to the pressure. Hence we omit the thermal expansion term, but retain the physical expansion term. 5. The tube wall is in equilibrium state at all times. The speed of propagation of dilatational waves in an elastic medium is denoted by s s E p 1 − νp λ p + 2µ p ve = ≡ ρp 1 + ν p 1 − 2ν p ρ p where ρ p is the polythene density, and λ p and µ p are Lame elastic constants. The value for ve is 1830 m/s. Hence the effect of the pressure change inside the tube is transmitted across the wall in a time that is short in comparison with the reaction time. Hence we can assume that the tube wall is in equilibrium state at all times. 6. The thermoelastic dissipation effect is negligible, such that the wall temperature can be determined independently of the deformation of the body. Consider the mechanical coupling term in the heat conduction equation 3λ p + 2µ p αT0 ∂εkk ∂T + = D p T,mm ∂t ρpC p ∂t where tensor notation has been used. By employing an order of magnitude argument, we see that the thermal dissipation effect is negligible if λ p [α] [εθθ ] 1. ρp C p Taking [εθθ ] ∼ [ p]/E p = 0.07, [C p ] = 2100 J/kg K, [ρ p ] = 920 kg/m3 , λ p = 2.96 × 109 Pa and α = 2.5 × 10−4 /K, we have λ p [α] [εθθ ] = 2.68 × 10−2 ρp C p which is small. Hence the thermal dissipation effect is negligible. 7. We assume that the tube does not degrade. This is a good assumption, because the reaction time scale is small; only a very thin region near the tube wall reaches the degradation temperature. 3. The governing equations Now we consider the one-dimensional propagation of a detonation wave down detonation transmission tubing of length L. Let the end of the tube at x = 0 be open and the end of the tube at x = L be closed. Let γi (x, t) denote the concentration of one of the phases in the tube where, here and with other variables, the subscripts 1, 2, 3 are used to indicate HMX, Al and inert gas, respectively. For solid materials on the tube wall, the (volume) concentration is taken as the average density across the tube. Let e(x, t) denote the specific internal energy, u(x, t) the gas speed, T (x, t) the absolute gas temperature, a(x, t) the inner radius of the tube, and p(x, t) the gas pressure. In this study, we assume that the tube radius is a linear function of pressure. Consider a small section of tube and let a0 and b0 be the initial inner and outer radii of the unstressed tube. We assume that the tube wall is in an equilibrium state with the inside and outside pressures and that the tube is displaced only radially at any location. Let pi and p0
722
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
be the uniform internal and external pressures, respectively. Then the general solution for the small displacement δ at r = a0 is ! # " pi a02 − p0 b02 a0 b02 ( p0 − pi ) 1 − νp 1 + νp δ= a0 − . (4) Ep Ep b02 − a02 b02 − a02 Thus we can write the inner radius of tube a as a = a0 + δ = Z 1 + Z 2 p,
(5)
where Z 1 = a0 −
2a0 b02 p0 E p b02 − a02
(6)
and Z2 =
Ep
h i a0 a02 1 − ν p + b02 1 + ν p . 2 2 b0 − a0
(7)
For the gas/particles motion inside the tube, the continuity equations for the aluminium can be written as ∂ 2 ∂ 2 a γ2 + a uγ2 = a 2 θΓ . ∂t ∂x
(8)
The first and second terms on the left-hand side are, respectively, the rate of change in, and the net mass flux of the aluminium into, the control volume. The right-hand side is a source term due to chemical reaction. The factors of a 2 are generated by the proportionality of the volume element to the cross-sectional area of the tube. Similarly, the continuity equation for the gas can be written as ∂ 2 ∂ 2 a γ3 + a uγ3 = a 2 (1 − θ) Γ . ∂t ∂x
(9)
In Eqs. (8) and (9), we have used the assumption that the gas and the aluminium flakes move with the same speed, i.e. u = u2 = u3. We assume that the gas/particles mixture is in pressure equilibrium and the overall momentum equation for the aluminium flake and the inert gas is i i ∂ h 2 2 ∂p ∂ h 2 a u (γ2 + γ3 ) + a u (γ2 + γ3 ) = −a 2 . ∂t ∂x ∂x
(10)
It is interesting to note that the momentum equation (10) cannot be written in conservation form due to the term on the right-hand side. It is not desirable to have an equation in non-conservation form. Now the right-hand side of Eq. (10) can be written as a2
∂p ∂a 2 p ∂a 2 ∂a 2 p ∂p = −p = − 2a Z 2 p ∂x ∂x ∂x ∂x ∂x
(11)
˜ x = L x, ˜ and p = [ p] p, ˜ where Eq. (5) has been used. If we non-dimensionalize Eq. (11) with the scalings a = a0 a, where [ p] = 12 MPa is the peak pressure, then the equation becomes ∂ p˜ ∂ a˜ 2 p˜ 2π p˜ ∂ p˜ a˜ 2 = − a˜ 2 (12) ∂ x˜ ∂ x˜ 1 + π p˜ ∂ x˜ where π = [ p]Z 2 /a0 equals 0.13. Thus, since π is still small even at the peak pressure, the right-hand side of Eq. (10) can be approximately written as a2
∂a 2 p ∂p ≈ ∂x ∂x
(13)
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
and the momentum equation becomes i i ∂ h 2 ∂ h 2 2 a u (γ2 + γ3 ) + a u (γ2 + γ3 ) + a 2 p = 0. ∂t ∂x The energy equation of the mixture is 2 2 u u ∂ a 2 γ1 e1 + a 2 γ2 + e2 + a 2 γ3 + e3 ∂t 2 2 2 2 ∂ a 2 pu ∂ u u + a 2 γ2 u =− + e2 + a 2 γ3 u + e3 − 2aΨ ∂x 2 2 ∂x
723
(14)
(15)
where Ψ is the heat loss rate per unit surface area. We now assume that the internal energies are linear functions of temperature with constant specific heat. We may express this as e1 = e10 + C1 (T − T0 ) , e2 = e20 + C2 (T − T0 ) , e3 =
e30
(16)
+ C3 (T − T0 ) ,
where e10 = θe20 + (1 − θ )e30 + 1H , e20 and e30 denote the respective internal energies at ambient room temperature, and −1H is the heat of reaction. There is no loss of generality in setting e20 = e30 = 0 and the Eq. (15) can be written as 2 2 u ∂ u 2 2 + e2 + a γ3 + e3 a γ2 ∂t 2 2 2 2 u u ∂ 2 2 2 + a γ2 u + e2 + a γ3 u + e3 + a up = a 2 Γ 1H − 2aΨ . (17) ∂x 2 2 The reaction equation is ∂ a 2 γ1 = −a 2 Γ ∂t where ( 2 ∗ 3 Γ = φγ1 T ≥ T ∗, 0 T
(18)
(19)
1
and φ = 1.7 × 105 kg 3 /ms is a constant. The derivation of φ is discussed in the Appendix. We assume that the gas obeys the perfect gas law p ≈ γ3 RT,
(20)
where R = Rgas /W , Rgas is the gas constant and W is the average molecular weight of the inert gas. We have made an assumption that the gas density ρ3 ≈ γ3 , since the solid phases occupy a very small volume. For easy presentation, it is more convenient to write the above flow equations in the vector form ∂v ∂F (v) + = S (v) ∂t ∂x
(21)
where a 2 γ2 a 2 γ3 v= (γ2 +γ3 ) a 2 u 1 2 1 2 2 a γ2 e2 + u + γ3 e3 + u 2 2
(22)
724
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
ua 2 γ2 ua 2 γ3 2 2 F (v) = a p (γ2 + γ3 ) a u + 1 2 1 2 2 a u γ2 e2 + u + γ3 e3 + u + p 2 2
(23)
a 2 θΓ a 2 (1 − θ) Γ . S (v) = 0 2 a Γ 1H − 2aΨ
(24)
and
Eq. (21) is a system of hyperbolic conservation laws and may be converted to an equivalent set of wave equations that represent nonlinear waves propagating at characteristic speeds, which in turn are functions of the local solutions, and vary in space and time. These characteristic speeds are λ 1 = u − cs ,
λ2 = λ3 = u,
λ4 = u + cs
(25)
where cs is the speed of sound: cs2 = B1 B2 , 1 , B1 = 2 a (γ2 + γ3 )
Rγ3 B2 = a 2 p 1 + C2 γ2 + C3 γ3
! .
(26)
It is also useful to express the equations in terms of the four primitive variables, γ2 , γ3 , u and p. This facilitates the conversion of the boundary conditions to a form that is more convenient for numerical computational purposes. In terms of these variables, the equations are ∂ a 2 γ2 = −d1 + Y1 , ∂t ∂ a 2 γ3 = −d2 + Y2 , ∂t ∂u = −d3 + Y3 , ∂t ∂ a2 p = −d4 + Y4 ∂t
(27)
where a 2 γ2 (Υ1 + Υ4 ) + Υ2 , B2 a 2 γ3 d2 = (Υ1 + Υ4 ) + Υ3 , B s2 B1 d3 = (Υ4 − Υ1 ) , B2 d1 =
d4 = Υ1 + Υ4 ,
(28)
725
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
# ∂ a2 p cs ∂u − , ∂x B1 ∂ x " # ∂ a 2 γ2 a 2 γ2 ∂ a 2 p Υ2 = λ2 − , ∂x B2 ∂x " # ∂ a 2 γ3 a 2 γ3 ∂ a 2 p − , Υ3 = λ3 ∂x B2 ∂x " # λ4 ∂ a 2 p cs ∂ (u) Υ4 = + , 2 ∂x B1 ∂ x λ1 Υ1 = 2
"
(29)
and Y1 = a 2 θΓ , Y2 = a 2 (1 − θ) Γ , uΓ , Y3 = − γ2 + γ3 Y4 =
Rγ3 γ2 C2 + γ3 C3
(30) !
p u2 a 2 Γ p (1 − θ) . a Γ 1H − [θ (C2 − C3 ) + C3 ] − T0 + − 2aΨ + 2 γ3 Rγ3 2
Finally, we denote the temperature within the tube wall as ψ. Inside the tube wall we have the heat diffusion equation in cylindrical coordinates (r, θ, x), 2 ∂ ψ ∂ψ ∂ψ 1 ∂ = Dp r (31) + ∂t r ∂r ∂r ∂x2 where D p = 1.03 × 10−7 m2 /s is the thermal diffusivity of the wall. Let y = r − a0 be the distance measured from the tube wall, where a0 = 1.5 mm denotes the tube radius. Now the length scale that the heat penetrates into the tube is p [y] = D p tr = 2.03 × 10−6 m where tr = 40 µs is the time scale for reaction. We can see that [y] a0 . This allows us to approximate the second term on the right-hand side of (31) by ∂ 2 ψ/∂ y 2 . Also, since the length scale for reaction lr = 0.08 m is much larger than [y], we can neglect the axial conduction term, ∂ 2 ψ/∂ x 2 . The heat diffusion equation (31) thus becomes ∂ψ ∂ 2ψ = Dp 2 , ∂t ∂y
(32)
where ψ = ψ(x, y, t). Note also that the length scale [y] of penetration is a lot less than the thickness h of the tube, since h = 8.5 × 10−4 m. Thus, the far region of the tube wall effectively remains unheated as the detonation wave passes. As a consequence of this, we are able later to reduce the domain of integration for the wall and hence speed up computations. 4. Numerical method 4.1. Initial and boundary conditions Each eigenvalue λi obtained in the previous section represents the characteristic speed at which a particular wave mode propagates. At a point on the boundary, some of the characteristics describe outgoing waves, while some of them describe incoming waves. The behaviour of the outgoing waves is completely determined by data contained
726
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
within and on the boundary of the computational domain. The number of boundary conditions that must be specified at a point on the boundary is equal to the number of incoming waves at that point. We specify the boundary conditions that determine the value of Υi for incoming waves, and compute, from the definition (29), the value of Υi for outgoing waves. At the open end, the gas is free to flow out of the tube. The number and type of boundary conditions required may vary from time to time at this boundary. There are two possible cases in this study. (1) Supersonic outflow: supersonic outflow at the open end is characterised by u > cs , so that all λi are greater than zero. Consequently, we can specify no boundary conditions at all, and the evolution of the flow at the boundary is determined completely by interior data. (2) Subsonic outflow: subsonic outflow is characterised by 0 > u > −cs , for which λ1−3 < 0 and λ4 > 0. Thus, we compute Υ1−3 from their definitions in (29) and specify Υ4 according to some boundary condition. In this study we specify the pressure to be one atmosphere. We achieve this by requiring that ∂(a 2 p)/∂t = 0. This condition is satisfied if Υ4 + Υ1 − Y4 = 0, giving Υ1 = Y4 − Υ4 . The other end of the tube is closed. The initial state has zero speed at the closed √ end, and the boundary conditions are chosen to keep the speed zero at all time. This condition is satisfied if − B1 /B2 (Υ4 − Υ1 ) + Y3 = 0, giving √ Υ1 = B2 /B1 Y3 + Υ4 . We choose the initial conditions to simulate hot combustion gases over a short distance L 1 near the free boundary, i.e. 0 ≤ x ≤ L 1 L. The gas pressure is taken to be constant throughout the tube at the external ambient value. A linear variation of gas temperature is assumed to be present near the open tube end: x , 0 ≤ x ≤ L 1, T0 + (Ti − T0 ) 1 − (33) T (x, 0) = L1 T0 , L1 < x ≤ L , where T0 is the reference temperature (nominally 293 K) and L 1 equals 0.1 m. This suffices to initiate the combustion process inside the tube, provided that the temperature Ti resulting from the ignition source exceeds the ignition criterion for the HMX. For the tube wall, we assume that there is no heat loss to the surrounding medium, so we impose an insulation condition at infinity. At the inner tube wall boundary, we need to specify two boundary conditions. These are continuity of the temperature field ψ (x, 0, t) = T (x, t) , and continuity of heat flux ∂ψ −Ψ = K p ∂ y y=0
(34)
(35)
where K p = D p C p ρ p is the wall conductivity, C p is the specific heat of the tube, and ρ p is density of the tube. As usual, the wall is assumed to be initially at ambient temperature throughout. This may be written as ψ (x, y, 0) = T0 .
(36)
4.2. Numerical solution procedure We use the time splitting method to obtain the numerical solution of Eq. (21). The first step is to remove the inhomogeneous term S(v) for Eq. (21). Thus the system reduces to ∂v ∂F (v) + = 0. ∂t ∂x
(37)
Eq. (37) can be solved using the TVD MacCormack method [16]. Let v¯ be the numerical solution of Eq. (37). Then the system of ordinary differential equations given by dv = S (¯v) , dt
(38)
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
727
together with the reaction Eq. (18), are solved in turn. Eq. (38) is solved using a Runge–Kutta Scheme by Cash and Krap [17], by utilizing the solution v¯ of Eq. (37) to determine the inhomogeneous term S. The heat conduction equation in the wall is parabolic in the time variable. We solve it by sequential application of the method of lines, integrating in any one step from time tn to tn+1 . We define 1y to be the mesh spacing within the wall. Let m 1 denote the number of x-mesh spacings and m 2 denote the number of y-mesh spacings. Then m 1 = L/1x and m 2 ≤ h/1y. We introduce mesh functions defined as the solution of the finite difference equations below and which are related to the exact solution by ψi, j ' ψ(i1x, j1y, t), where 0 ≤ i ≤ m 1 , 0 ≤ j ≤ m 2 . Thus, using standard three-point central differencing, we obtain the following ordinary differential equation (ODE): ψi, j+1 − 2ψi, j + ψi, j−1 dψi, j = . dt (1y)2
(39)
At each time step, we set ψi,0 = T (i1x, tn ) and ψi,m 2 = T0 . Since the heat is only absorbed in a very thin region near the tube wall, we choose m 2 to equal one tenth of h/1y, i.e. we only compute the temperature over one tenth of the tube thickness in our computation. This saves a lot of computational time. Eq. (39) consists of a system of ODEs that are solved using the same Runge–Kutta method. Thus, the inhomogeneous term and the heat loss term are solved together. 5. Numerical results To demonstrate the importance of the tube wall interactions, three sets of results will be presented here. The first set is for results without any tube wall interactions, the second set of results is with wall heat loss only, and the third set is for the complete model. All the solutions are depicted at equal time intervals of 0.1 ms. The tube length is L = 2 m and 801 nodes are used in the calculation. 5.1. Without any tube wall interactions In this section, the numerical results without any tube wall interactions are presented, i.e. the radius a is a constant and Ψ = 0. The graphs of γ1 , γ2 , γ3 , p, T and Mach number as a function of x are shown in Fig. 1. The very last profile shows conditions after the shock has hit the closed end of the tube and reflected. The HMX concentration Fig. 1(a) shows that the reaction zone extends about 0.11 m. The aluminium flake and gas concentration are shown in Fig. 1(b) and (c), respectively. We can see that they have similar features. This is because of our assumption that the gas and aluminium flake move with the same speed. The peak pressure in the pressure Fig. 1(d) is about 40 MPa and the pressure at the open end stays nearly constant over the calculation. Fig. 1(e) depicts the solution for the temperature. The peak temperature is about 5300 K. The Mach number Fig. 1(f) shows that the gas exits at sonic speed and the Mach number near the front is about 0.7. The detonation speed is about 2800 m/s. From our numerical results, the peak pressure is about 40 MPa and the detonation speed is about 2800 m/s. However, if we compare these with information from ICI, we find that our peak pressure is about four times more than the ICI observations. Also, our numerical detonation speed is faster than their experimental observation of 2000 m/s. 5.2. The effect of heat loss to the tube wall In this subsection, we have included the heat loss term in our model. The numerical solutions are shown in Fig. 2. The pressure result in Fig. 2(a) shows a sharp discontinuity. The peak pressure is about 18 MPa, which is a reduction of over 50% in magnitude compared to the result without heat loss. The detonation speed is about 2000 m/s, which is now in good agreement with the ICI experiment results. The gas temperature Fig. 2(b) shows that the peak temperature is about 3200 K, which is about half of the maximum temperature rise without the heat loss effect. The aluminium flake and gas concentration are shown if Fig. 2(c) and (d), respectively. We can see that they have smaller peak values in comparison with the results without heat loss. Fig. 3 shows the temperature inside the wall. As we expected, the heat is only absorbed within a very thin region near the interior tube wall.
728
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
(a) HMX concentration.
(b) Al/Al2 O3 concentration.
(c) Gas concentration.
(d) Pressure.
(e) Gas temperature.
(f) Mach number. Fig. 1. Numerical results without any tube wall interactions at equal time intervals of 0.1 ms.
5.3. Full model results The full model results are shown here. The graph of pressure is shown in Fig. 4(a). We see that the peak pressure is about 13 MPa. This is about a one third reduction compared with the previous result. Physically, this is because the gas now occupies a bigger volume. The corresponding temperature field is shown in Fig. 4(b). There is not much difference between this and the temperature graph without tube expansion from the previous section. This is because we have no energy loss from the tube expansion. The graph of the inner tube radius is shown in Fig. 4(c). We see that it shows a similar profile to the pressure graph. This is because they are connected by the linear relation equation (5). There is about a 14% increase in the tube radius at the peak pressure.
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
(a) Pressure.
(c) Al/Al2 O3 concentration.
729
(b) Gas temperature.
(d) Gas concentration. Fig. 2. Numerical results for heat loss model at equal time intervals of 0.1 ms.
Fig. 3. Tube wall temperature profile at 0.8 ms.
The solutions for the gas concentration and aluminium particles concentration are given in Fig. 4(d) and (e), respectively. We see that, as before, there is a constant state in between the detonation and the rarefaction. Also,
730
(a) Pressure.
(c) Inner tube radius.
(e) Al/Al2 O3 concentration.
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
(b) Gas temperature.
(d) Gas concentration.
(f) Gas speed. Fig. 4. Numerical results for full model at equal time intervals of 0.1 ms.
it is more noticeable than before. The gas speed is shown in Fig. 4(f). The typical gas speed is about 800 m/s after the shock wave. 6. Discussions Our model has a ZND wave structure, that is, a shock followed by a deflagration wave in the HMX, including the tube wall heat loss effect and the tube expansion effect. The numerical results of our model compare favourably with experiments performed by ICI. However, there are some aspects of the results that we need to discuss further. If the pressure results in Fig. 4(a) are examined carefully, a pulsation detonation pressure profile is observed. Fig. 5 shows the maximum detonation pressure against time, which clearly reveals that the model is unstable. The frequency
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
731
Fig. 5. Maximum detonation pressure against time (dashed line, 801-nodes; solid line 1601-nodes).
of the oscillation remains approximately constant when we enlarge the number of nodes. Thus we believe that the instability is not caused by the methodology. As we mentioned at the beginning of this paper, the ZND wave structure is usually unstable [13] and our numerical calculations clearly show a one-dimensional instability. Another interesting feature is the constant state in between the reaction zone and the head of the rarefaction wave: in the first set of results the double heads emerge in the solution. This is due to the final state of the detonation moving too fast such that it runs away from the following rarefaction wave. We can see from the pressure graph Fig. 1(d) that the first head begins to develop at 0.4 ms and its strength begins to increase, whereas the second head strength starts to decrease. Also, the distance between the two heads increases and the detonation moves further away from the rarefaction wave. The double-head feature disappeared in the second set of results, but there is a distinguishable region in between the detonation and the head of the rarefaction wave. The wall heat loss effect caused some energy in the gas to be transferred to the tube wall. This reduces the strength of the detonation front and the rarefaction wave. However, heat loss has much effect on the rarefaction wave and reduces the wave propagation speed relative to the detonation. Hence the second head vanishes in the solution and we can see a constant state. This constant state becomes more noticeable when we add the tube expansion effect. The speed of propagation of the rarefaction waves in between the second and third set of results is about the same. But this time the tube expansion significantly reduces the detonation head strength, such that we see the constant state more clearly, although there are some wrinkles along the constant state. We believe that this is caused by the instability of the model. This constant state tells us that our detonation speed is faster than the CJ speed and our final result has a weak solution. 7. Conclusions In this paper we have formulated a one-dimensional model in which the detonation wave is propagated in detonation transmission tubing and where the reaction rate is described by a single-step kinetic equation. We then went on to describe the characteristic and the numerical method of the model equation. The heat loss effect and the tube expansion effect have also been considered. The final results obtained are in good agreement with the experimental data. The results show that the heat loss and the tube expansion have a major effect on the structure of the detonation wave. The long reaction zone results in a non-adiabatic combustion system. The elastic material of tubing expands under high interior pressure. Also, the constant state in between the detonation and the rarefaction wave tell us that we have a weak detonation. An instability feature causes an unstable (steady oscillation) detonation. Acknowledgements The author would like to express his gratitude to ICI for their generous financial support and to Dr A. Janes for his valuable help during this project. The author would also like to acknowledge the highly useful discussions about this problem with Dr. D. Sutton and Dr. P. Lynch of ICI Explosives. Appendix In this appendix, we derive and discuss the various physical scalings that are quoted in the text and used to justify the various assumptions made.
732
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
One physical quantity needed is the time scale of reaction. We were provided with the results of some experiments on DTT performed by ICI Explosives Division which enabled us to make an estimate of this. In one of these experiments, the light emitted by the detonation was measured as it passed a fixed point of the tube. It seems natural to assume that this light is generated by the reaction zone, and that the light zone τ is the required reaction time. The reaction constant φ in Eq. (18) can then be chosen to generate such a reaction by choosing φ = 3(γ10 )1/3 /τ , where γ10 = 11.09 kg/m3 is the initial concentration of HMX. The observed light period is about τ = 100 µs. The use of this value in our numerical calculations, however, predicted a detonation speed significantly below the experimentally observed speed of 2 km/s. Also, the intensity of the light emission profile does not fit the theoretical profile well in the early stages. There is a 15 µs lag before the maximum light output was achieved, whereas the theoretical profile has its maximum attained immediately. Such disagreements are not entirely unexpected, given the highly simplified assumptions made about the kinetics and the geometry. Eventually, we elected to choose a value for τ that gave good agreement with the observed detonation speed, rather than with the light emission, and the value satisfying this criterion, τ = 40 µs, was found by numerical experiments. One of the assumptions made was that the speed of a typical aluminium particle was the same as that of the gas flow. A typical Reynolds number for an aluminium flake is ρulal = 1691 µ
Ral =
where u ∼ 800 m/s is a typical gas speed (taken from the final numerical results), ρ ∼ 16 kg/m3 is a typical gas density after detonation (also from the numerical results), µ ∼ 7.57 × 10−5 kg/ms is the dynamic gas viscosity (of air) at 2500 K, and lal = 10 µm is the relevant length scale for gas flow past a thin plate. Since the Reynolds number Ral is large, the motion of the aluminium obeys inviscid dynamics. In inviscid flows the acceleration of the particle is caused by the pressure difference on the forward and backward faces of the particle. This pressure difference has the order of magnitude 21 ρu 2 , so let us take a mean pressure difference as 1 1 2 1p = ρu . 2 2 For an aluminium flake, this acts over an area of lal2 /2 (because of the aspect ratio effect). Therefore, the typical force is Fal =
lal2 ρu 2 . 8
By Newton’s second law, Mal
du al = Fal , dt
where Mal is the mass of each aluminium particle. So the time 1tal needed to achieve a speed change 1u al = u is 1tal =
Mal u 8ρal dal = = 1.68 × 10−7 s, Fal ρu
where dal = 0.1 µm is the thickness of the aluminium flake and ρal is the density of aluminium. We can see that the time scale is fast compared with the reaction time scale tr = 40 µs. Hence we assume that aluminium particles move with the speed of the gas. A similar calculation for HMX, which has density ρHMX = 1900 kg/m3 , gives 1tHMX = 2.71 × 10−5 s. We have assumed that the HMX is stationary, but clearly this not as good an approximation as the last one. The particles would be accelerated to the gas speed over a time comparable to the reaction time, and ideally the model should be modified. However, the chemical reaction is proportional to the surface area, so most of the reaction will take place during the early stages of the acceleration process when the speed that has been achieved is small. Hence we expect the effect of the motion to be relatively minor, and so neglect it.
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
733
The time scale to heat up a HMX particle by diffusion is 2 1 2 lHMX = 4.0 × 10−3 s, [t] = DHMX where DHMX = 1.0 × 10−7 m2 /s is the diffusivity of the HMX. This is slow compared to the time scale for the reaction. We can conclude from this that the bulk of the particle does not have time to heat up by diffusion, and hence that it does not react as a uniform whole, and must therefore react by conversion of the surface of HMX to gas. Heat diffusion takes place, but only in a thin layer next to the newly exposed surface. This is rapid because of the large temperature gradient. Diffusion into the bulk is negligible. We thus assume in the model that the rate of reaction Γ is proportional to the surface area of HMX remaining, and that the temperature of the unburnt HMX remains constant during the reaction. Conversely, the time scale taken to heat up the aluminium flake of thickness dal is 2 1 2 lal = 2.56 × 10−11 s, [t] = Dal where Dal = 9.76 × 10−5 m2 /s is the diffusivity of the aluminium. Now this value is very short compared with our reaction time scale, so we assume that the temperature of the aluminium flake is the same as that of the gas. From thermodynamic principles, the oxygen reacts with the metal preferentially to form metal oxide, then any remaining oxygen reacts with carbon to form carbon monoxide, and then the remnant oxygen reacts with hydrogen to form water. According to this rule, we can write the overall chemical formula per unit mole of HMX as Eq. (1). We can see from Eq. (1) that the reaction product is dominated by carbon monoxide and molecular nitrogen. We treat the gases as a single gas product, taking the molecular weight of the hypothetical gas to be the mole average of the gas products. The average molecular weight W = 22.99 g/mol is used here. Since the aluminium flake is heated instantaneously, we assume that the oxygen initially present in the air reacts immediately with aluminium to form Al2 O3 , such that no free oxygen remains in the air after the passage of the shock. The chemical formula for this reaction is 0.954Al + 0.155O2 → 0.103Al2 O3 + 0.747Al. We see that only about 22% of that aluminium becomes Al2 O3 from this reaction. To avoid working with a two-reaction kinetic scheme, we decided to omit this reaction. For all the aluminium to become Al2 O3 , we need a further 0.56 mole of oxygen from the 1 mole of HMX. The chemical formula is 0.747Al + 0.560O2 → 0.375Al2 O3 . weight of oxygen So the mass ratio θ = 0.56×molecular = 0.061. We can see that only a very small fraction of gas is molecular weight of HMX absorbed by aluminium. Most of the HMX becomes gas.
References [1] M. Berthelot, P. Vielle, On the velocity of propagation of explosive processes in gases, C. R. Hebd. Sceances Acad. Sci. 93 (1881) 18–21. [2] M. Berthelot, P. Vielle, On explosive waves, C. R. Hebd. Sceances Acad. Sci. 94 (1882) 882. [3] E. Mallard, H. Le Chatelier, On the propagation velocity of burning in gaseous explosive mixtures, C. R. Hebd. Sceances Acad. Sci. 93 (1881) 145–148. [4] D.L. Chapman, On the rate of explosive in gases, Philos. Mag. 47 (1899) 90–104. [5] E. Jouguet, On the propagation of chemical reactions in gases, J. de Math. Pures et Appl. 1 (1905) 347–425. Continued in 2 (1906) 5 – 85. [6] E. Jouguet, Mecanique des Explosifs, Octave Doin et Fils, Paris, 1917. [7] Y.B. Zeldovich, On the theory of the propagation of detonation in gaseous system, Zh. Eksp. Teor. Fiz. 10 (1940) 542–568. English translation: NACA TM1261, 1960. [8] J. von Neumann, Theory of detonation wave, Technical Report, OSRD Report No. 549, Paris, 1917. [9] W. Doering, On detonation processes in gases, Ann. Phys. 43 (1943) 421–436. [10] A.K. Oppenheim, J. Rosciszewski, Determination of the detonation wave structure, in: Proc. 9th Symp. Int. on Combustion, 1963, pp. 424–441.
734
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
[11] D.H. Edwards, G.T. Williams, J.C. Breeze, Pressure and velocity measurements on detonation waves in hydrogen–oxygen mixtures, J. Fluid Mech. 6 (1959) 497–517. [12] J.H.S. Lee, Dynamic parameters of gaseous detonation, Ann. Rev. Fluid Mech. 16 (1984) 311–336. [13] G.E. Abouseif, T.Y. Toong, Theory of unstable one-dimensional detonations, Combust. Flame 45 (1982) 67–94. [14] F. Zhang, H. Gronig, Spin detonation in reactive particles-oxidizing gas flow, Phys. Fluids A 3 (1991) 1983–1990. [15] T.B. Brill, H. Arisawa, P.J. Brush, P.E. Gongwer, G.K. Williams, Surface chemistry of burning explosives and propellants, J. Phys. Chem. 99 (1995) 1384–1392. [16] S.F. Davis, TVD finite difference schemes and artificial viscosity, Technical Report, NASA CR 172373, 1984. [17] J.R. Cash, A.H. Karp, A variable order Runge–Kutta Method for initial value problems with rapidly varying right-hand sides, ACM Trans. Math. Software, 16 201–222.
A numerical study of a detonation wave in detonation transmission tubing D.K.L. Tsang ∗ School of Mechanical, Aerospace and Civil Engineering, The University of Manchester, PO Box 88, Sackville Street, Manchester M60 1QD, UK Received 12 July 2004; received in revised form 15 December 2005; accepted 9 February 2006
Abstract In this paper we formulate a one-dimensional detonation wave model for detonation transmission tubing. The reaction is described by a temperature-independent single-step kinetic equation and the resulting equations are a system of hyperbolic conservation laws. The initial boundary value problem is then solved by using the TVD MacCormack method. The effects of heat loss to the tube wall and of the tube expansions have also been considered. The result obtained shows that the heat loss and the tube expansion have a major effect on the pressure prediction. c 2006 Elsevier Ltd. All rights reserved.
Keywords: Detonation; Rarefaction wave; Detonation transmission tubing; Deflagration; HMX
1. Introduction Detonation transmission tubing (DTT) is a modern method for transmitting the detonation signal from a start device to the detonator delay element and hence to the main explosive. It is much safer than electric wiring, since conventional wire can be accidentally initiated by an electric discharge. The DTT consists of a tube formed from SURLYN, a copolymer of ethylene and methacrylic acid. The internal and external tube diameters are 1.3 mm and 3.0 mm, respectively. The interior surface is coated with a mixture of aluminium (Al) flakes and HMX. Core loadings are typically 16 mg/m (milligrams per meter length of the tube). About 92% of the mixture is the HMX and this is the main provider of energy to support the detonation wave. The rest is aluminium, which is mainly present to expedite initiation. Each HMX particle is approximately in the form of a cube with sides of length 40 micron. Each Al flake is in the form of a thin plate approximately 10 micron square and 0.1 micron thick. The mixture occupies about 1% of the internal tube volume. The DTT is normally initiated by a shock wave at one end. Associated with the shock wave are pressure, temperature and density changes. It is posited that the passage of the shock wave raises the temperature of the gas in the tube, and this subsequently causes the HMX to react. The temperature rise also heats the Al particles sufficiently for the oxide layer on the surface to soften and for rapid reaction between the exposed metal and oxygen to commence. There is interest at the present moment in trying to predict the effects of any accidental damage that may occur to the ∗ Tel.: +44 1612754506; fax: +44 1612754901.
E-mail address: [email protected]. c 2006 Elsevier Ltd. All rights reserved. 0895-7177/$ - see front matter doi:10.1016/j.mcm.2006.02.008
718
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
Nomenclature W ρal ρHMX ρp Dal DHMX Dp Kp C1 C2 C3 Cp 1H Rgas R L L1 T∗ T0 Ti p0 a0 b0 Ep νp
average molecular weight 22.99 g/mol density of aluminium 2700 kg/m3 density of HMX 1900 kg/m3 density of tube wall 920 kg/m3 diffusivity of aluminium 9.76 × 10−5 m2 /s diffusivity of HMX 1.00 × 10−7 m2 /s diffusivity of tube wall 1.03 × 10−7 m2 /s conductivity of tube wall 0.198 W/m K specific heat of HMX 1500 J/kg K specific heat of Al 903 J/kg K specific heat of gas 1000 J/kg K specific heat of tube wall 2100 J/kg K heat of reaction 5.66 × 106 J/kg gas constant 8.31 J/mol K perfect gas constant 361.65 J/kg K length of tube 2 m initial length of temperature variation 0.1 m cut-off temperature of HMX 450 K reference temperature 293 K ignition temperature 500 K reference pressure 1.01 × 105 Pa initial inner radius 6.5 × 10−4 m initial outer radius 1.5 × 10−3 m Young’s modulus of tube 1.80 × 108 Pa Poisson’s ratio of tube 0.49
tube during mining operations. For this purpose it would be useful to have a model of the operation of the tube and to understand the mechanism of the detonation process. Detonation wave is a combustion process initiated by a shock. Detonations were first investigated scientifically by Berthelot and Vielle [1,2], and by Mallard and Le Chatelier [3]. Hypotheses leading to mathematical predictions of the speed of propagation of the waves were given by Chapman [4] and Jouguet [5,6]. A good approximation of the usual detonation structure is a shock wave followed by a deflagration. This shock–deflagration structure of detonation was first proposed by Zeldovich [7], von Neumann [8] and Doering [9] and is referred to as the ZND wave structure. On the basic of numerical calculations of steady, planar detonation structure (for example, Oppenheim and Rosciszewski [10]) and of good experimental measurements performed mainly in the 1960s (for example, Edwards et al. [11]), it was widely believed that the ZND structure was representative of most real detonations. However, more recent studies established that this structure is usually unstable. In fact, these detonations generally possess a complex three-dimensional substructure involving transverse waves that produce a steady, planar wave only on average (see Lee [12]). In Abouseif and Toong [13], the mechanism of instability is identified. In their analysis, they demonstrated that the interaction between the irreversible temperature fluctuations and the finite reaction zone induces an oscillating energy source field, which then leads to shock perturbations and thereby detonation pressure oscillations. Detonation propagation in premixed substances is one of the basic aspects of combustion. The complexity of detonations is greatly increased if the pre-mixed system contains finely pulverized organic or metallic solids suspended in air or oxygen. Self-sustained detonations can be obtained in heterogeneous dust fuel-gas oxidizer mixtures, and the propagation mechanism is dominated by a transverse wave structure similar to those in homogeneous gas mixtures. These facts have been firmly established by the observations of Zhang and Gronig [14]. In their reactive corn starch particles/oxidizing gas system, they find that the momentum and temperature relaxations are dominant immediately
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
719
behind the shock front. There, the gas velocity and temperature jump at once to the new state corresponding to the shock strength, while the particle velocity and temperature increase continuously from the initial state. This causes the induction zone length to increase. Heterogeneous effects like pyrolysis, surface oxidation, etc, also cause the induction zone length to increase. Experimental observations for the DTT, provided to us by ICI Explosives Division, proffer that the detonation wave passes down the DTT with a typical speed of 2 km/s and a peak pressure of about 10–12 MPa. In the present work, we begin with the simplest hypothesis and formulate a one-dimensional wave model for the passage of the detonation wave along an unbroken length of DTT. The reactants are considered to be homogeneously distributed within the tube with the reaction rate described by a temperature-independent single-step kinetic equation. This results in a set of hyperbolic conservation equations that we solve numerically. The effect of heat loss to the tube wall and the effect of tube expansion have also been considered. The model showed good agreement with experimental observations. 2. Assumptions In our wave model, the tube wall interaction terms have been considered in the derivation. We now list and discuss the key assumptions made in the gas/particulars motion and in the tube wall of our one-dimensional model. The various time scales used in the discussion are derived in the Appendix. 2.1. Assumptions for gas/particulars motion The assumptions for the gas/particulars motion are as follows: 1. The gas flow is one-dimensional, so that all variables depend only on x, the distance down the transmission tube and on t, the time. In practice, when a detonation wave propagates in a tube, the front structure is intrinsically two- or threedimensional and non-steady. Also, the wall cooling effect, the wall expansion effect and the heterogeneous nature of the HMX prevents the problem from being one-dimensional. However the one-dimensional detonation wave is the simplest representative system, is significantly easier to solve, gives insight into the basic physics of the problem, and needs to be understood before investigations of any non-planar phenomena can properly begin. Therefore, we begin by considering the one-dimensional calculation. 2. The speed of the aluminium flake is the same as the gas speed, but the speed of the HMX is negligible, i.e. the HMX remains stationary at its original position. A typical relaxation time for the speed of the aluminium flake is 0.193 µs. We can see that this time scale is short compared to a reaction time scale of 40 µs. We can thus assume that the aluminium flake is instantaneously accelerated by the gas flow, and so always has the same speed as the gas itself. However, a typical speed relaxation time for the HMX particles is 27.1 µs. Clearly, the assumption that the HMX is stationary is not perfect. The particles would be accelerated to gas speed over a time comparable to the time of the chemical reaction. However, during the time that the conversion of HMX to gas is a maximum, the particle motion is small, so we would expect the HMX motion effect to be small. 3. The temperature of the HMX remains constant during the reaction, but the temperature of the aluminium flake is assumed to be the same as the gas temperature. The time scale to heat up an HMX particle by diffusion is 4.0 × 10−3 s. This is slow compared to the time scale for the reaction. We can conclude from this that the bulk of the particle does not have time to heat up by diffusion. Hence it does not react as a uniform whole, and must therefore react by conversion of the surface of HMX to gas. Consequently, the temperature of the unburnt HMX remains constant until reaction occurs. Conversely, the time scale taken to heat up an aluminium flake is 2.56×10−11 s. Now this value is very short compared with our reaction time scale, so we assume that the temperature of the aluminium flake is the same as that of the gas. 4. The rate of reaction Γ (defined as the rate of change of mass of HMX per unit volume of tube) is zero below a certain cut-off gas temperature T ∗ . Above that gas temperature, Γ is proportional to the surface area of HMX remaining. The kinetics for a slab of uniformly heated HMX has been investigated by [15], but these cannot be applicable here. We have already given reasons why the reaction is expected to be a surface reaction, rather than a volume
720
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
reaction, and for this reason we assume that the rate of reaction Γ is proportional to the surface area of HMX remaining. It is known that the HMX decomposes at 450 K, and we choose this as our cut-off temperature T ∗ . 5. Only a single chemical reaction takes place. The heat of this reaction is the combined heat of reaction of HMX decomposition and aluminium oxidization. The aluminium is treated as an inert phase which accumulates mass as the reaction proceeds. This assumption was made to reduce the number of chemical reactions that the eventual computer program would have to keep track of, and hence to produce a more tractable model. A number of simplifications follow from the assumption. The final state of chemical equilibrium can be established from thermodynamic considerations, and in the model it is assumed that the reaction progresses evenly towards this state in proportion to the amount of HMX that has decomposed. In the real process, oxygen is likely to combine preferentially with the aluminium present. Only when this has been completely converted to Al2 O3 would there be significant carbon monoxide and water products. In the same way, the oxygen initially present in the tube would probably combine very early on, rather that continuously over the lifetime of the HMX decomposition. Thus, the real chemical reaction would be multi-stage. 6. During reaction, 1 unit mass of HMX produces θ unit masses of oxygen and 1 − θ unit masses of inert gas, where θ is a constant. The θ unit masses of oxygen then react with the aluminium instantaneously, such that no free oxygen is present in the model. The chemical kinetics of HMX/Al system are very complex and unknown to us. In order to produce a tractable model, we assume only a temperature-independent single-step reaction and that the heat of this reaction is the combined heat of HMX decomposition and Al oxidization. We write the overall chemical formula per unit mole of HMX as follows: HMX + 0.954Al + 0.704N2 + 0.155O2 → 0.477Al2 O3 + 4CO + 2.879H2 O + 1.121H2 + 4.704N2 .
(1)
We treat the gases as a single gas product, taking the W molecular weight of the hypothetical gas to be the mole average of the final gas products, equal to 22.99 g/mol. The mass ratio θ is equal to 0.061. We can see that only a very small fraction of gas is absorbed by the aluminium flake. Most of the HMX becomes gas. 2.2. Assumptions for tube wall The assumptions for the tube wall are as follows: 1. The tube deformation is small. This assumption implies that the displacements are sufficiently small so that no distinction is needed between the coordinates of a particle before and after deformation, and that the displacement gradients are sufficiently small so that their products may be neglected. 2. The material behaves elastically at all times. This assumption implies that neither the temperature changes nor the stresses are too large. However, in our present study there is a huge difference between the internal pressure and external pressure. From our numerical results, we find that there is about a 14% increase in the inner tube radius at the peak pressure. So the stress changes are large and the deformation is not small. Nevertheless, in order to produce a tractable model, we have made this assumption. 3. The tube is displaced only radially at any locations. This is a plane strain assumption, and applies when the length of the body is large compared with the maximum cross-sectional dimension. In this study, we have a long tube with very small cross-sectional area, so the concept of plane strain may be used. 4. The thermal expansion is negligible. If we consider the plane strain problem, then one of the stress–strain relations and one of the strain–displacement relations are 1 e e εθθ = σθθ − νe σrr + αT (2) Ee ue εθθ = (3) r
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
721
where ε is the strain, σ e is the stress, u e is the radial deformation displacement, r is the radial coordinate, E p = 1.8 × 108 Pa is the Young’s modulus, ν p = 0.49 is the Poisson ratio, and α = 2.5 × 10−4 /K is the thermal expansion coefficient. From an order of magnitude analysis, we have [ p] + α [T ] = 0.07 + 0.75 [εθθ ] ∼ Ep where [ p] = 12 MPa is the peak pressure and [T ] = 3000 K is the peak temperature. We can see that the strain is dominated by the thermal expansion and that, at r = a0 , the deformation displacement [u] = 0.82a0 by using Eq. (3). However, such a large displacement is not observed. Presumably this is because the heated layer is above its melting point and undergoes molecular rearrangement, rather than elastic displacement. If this is correct, we would expect to see no thermal expansion at all but only physical expansion due to the pressure. Hence we omit the thermal expansion term, but retain the physical expansion term. 5. The tube wall is in equilibrium state at all times. The speed of propagation of dilatational waves in an elastic medium is denoted by s s E p 1 − νp λ p + 2µ p ve = ≡ ρp 1 + ν p 1 − 2ν p ρ p where ρ p is the polythene density, and λ p and µ p are Lame elastic constants. The value for ve is 1830 m/s. Hence the effect of the pressure change inside the tube is transmitted across the wall in a time that is short in comparison with the reaction time. Hence we can assume that the tube wall is in equilibrium state at all times. 6. The thermoelastic dissipation effect is negligible, such that the wall temperature can be determined independently of the deformation of the body. Consider the mechanical coupling term in the heat conduction equation 3λ p + 2µ p αT0 ∂εkk ∂T + = D p T,mm ∂t ρpC p ∂t where tensor notation has been used. By employing an order of magnitude argument, we see that the thermal dissipation effect is negligible if λ p [α] [εθθ ] 1. ρp C p Taking [εθθ ] ∼ [ p]/E p = 0.07, [C p ] = 2100 J/kg K, [ρ p ] = 920 kg/m3 , λ p = 2.96 × 109 Pa and α = 2.5 × 10−4 /K, we have λ p [α] [εθθ ] = 2.68 × 10−2 ρp C p which is small. Hence the thermal dissipation effect is negligible. 7. We assume that the tube does not degrade. This is a good assumption, because the reaction time scale is small; only a very thin region near the tube wall reaches the degradation temperature. 3. The governing equations Now we consider the one-dimensional propagation of a detonation wave down detonation transmission tubing of length L. Let the end of the tube at x = 0 be open and the end of the tube at x = L be closed. Let γi (x, t) denote the concentration of one of the phases in the tube where, here and with other variables, the subscripts 1, 2, 3 are used to indicate HMX, Al and inert gas, respectively. For solid materials on the tube wall, the (volume) concentration is taken as the average density across the tube. Let e(x, t) denote the specific internal energy, u(x, t) the gas speed, T (x, t) the absolute gas temperature, a(x, t) the inner radius of the tube, and p(x, t) the gas pressure. In this study, we assume that the tube radius is a linear function of pressure. Consider a small section of tube and let a0 and b0 be the initial inner and outer radii of the unstressed tube. We assume that the tube wall is in an equilibrium state with the inside and outside pressures and that the tube is displaced only radially at any location. Let pi and p0
722
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
be the uniform internal and external pressures, respectively. Then the general solution for the small displacement δ at r = a0 is ! # " pi a02 − p0 b02 a0 b02 ( p0 − pi ) 1 − νp 1 + νp δ= a0 − . (4) Ep Ep b02 − a02 b02 − a02 Thus we can write the inner radius of tube a as a = a0 + δ = Z 1 + Z 2 p,
(5)
where Z 1 = a0 −
2a0 b02 p0 E p b02 − a02
(6)
and Z2 =
Ep
h i a0 a02 1 − ν p + b02 1 + ν p . 2 2 b0 − a0
(7)
For the gas/particles motion inside the tube, the continuity equations for the aluminium can be written as ∂ 2 ∂ 2 a γ2 + a uγ2 = a 2 θΓ . ∂t ∂x
(8)
The first and second terms on the left-hand side are, respectively, the rate of change in, and the net mass flux of the aluminium into, the control volume. The right-hand side is a source term due to chemical reaction. The factors of a 2 are generated by the proportionality of the volume element to the cross-sectional area of the tube. Similarly, the continuity equation for the gas can be written as ∂ 2 ∂ 2 a γ3 + a uγ3 = a 2 (1 − θ) Γ . ∂t ∂x
(9)
In Eqs. (8) and (9), we have used the assumption that the gas and the aluminium flakes move with the same speed, i.e. u = u2 = u3. We assume that the gas/particles mixture is in pressure equilibrium and the overall momentum equation for the aluminium flake and the inert gas is i i ∂ h 2 2 ∂p ∂ h 2 a u (γ2 + γ3 ) + a u (γ2 + γ3 ) = −a 2 . ∂t ∂x ∂x
(10)
It is interesting to note that the momentum equation (10) cannot be written in conservation form due to the term on the right-hand side. It is not desirable to have an equation in non-conservation form. Now the right-hand side of Eq. (10) can be written as a2
∂p ∂a 2 p ∂a 2 ∂a 2 p ∂p = −p = − 2a Z 2 p ∂x ∂x ∂x ∂x ∂x
(11)
˜ x = L x, ˜ and p = [ p] p, ˜ where Eq. (5) has been used. If we non-dimensionalize Eq. (11) with the scalings a = a0 a, where [ p] = 12 MPa is the peak pressure, then the equation becomes ∂ p˜ ∂ a˜ 2 p˜ 2π p˜ ∂ p˜ a˜ 2 = − a˜ 2 (12) ∂ x˜ ∂ x˜ 1 + π p˜ ∂ x˜ where π = [ p]Z 2 /a0 equals 0.13. Thus, since π is still small even at the peak pressure, the right-hand side of Eq. (10) can be approximately written as a2
∂a 2 p ∂p ≈ ∂x ∂x
(13)
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
and the momentum equation becomes i i ∂ h 2 ∂ h 2 2 a u (γ2 + γ3 ) + a u (γ2 + γ3 ) + a 2 p = 0. ∂t ∂x The energy equation of the mixture is 2 2 u u ∂ a 2 γ1 e1 + a 2 γ2 + e2 + a 2 γ3 + e3 ∂t 2 2 2 2 ∂ a 2 pu ∂ u u + a 2 γ2 u =− + e2 + a 2 γ3 u + e3 − 2aΨ ∂x 2 2 ∂x
723
(14)
(15)
where Ψ is the heat loss rate per unit surface area. We now assume that the internal energies are linear functions of temperature with constant specific heat. We may express this as e1 = e10 + C1 (T − T0 ) , e2 = e20 + C2 (T − T0 ) , e3 =
e30
(16)
+ C3 (T − T0 ) ,
where e10 = θe20 + (1 − θ )e30 + 1H , e20 and e30 denote the respective internal energies at ambient room temperature, and −1H is the heat of reaction. There is no loss of generality in setting e20 = e30 = 0 and the Eq. (15) can be written as 2 2 u ∂ u 2 2 + e2 + a γ3 + e3 a γ2 ∂t 2 2 2 2 u u ∂ 2 2 2 + a γ2 u + e2 + a γ3 u + e3 + a up = a 2 Γ 1H − 2aΨ . (17) ∂x 2 2 The reaction equation is ∂ a 2 γ1 = −a 2 Γ ∂t where ( 2 ∗ 3 Γ = φγ1 T ≥ T ∗, 0 T
(18)
(19)
1
and φ = 1.7 × 105 kg 3 /ms is a constant. The derivation of φ is discussed in the Appendix. We assume that the gas obeys the perfect gas law p ≈ γ3 RT,
(20)
where R = Rgas /W , Rgas is the gas constant and W is the average molecular weight of the inert gas. We have made an assumption that the gas density ρ3 ≈ γ3 , since the solid phases occupy a very small volume. For easy presentation, it is more convenient to write the above flow equations in the vector form ∂v ∂F (v) + = S (v) ∂t ∂x
(21)
where a 2 γ2 a 2 γ3 v= (γ2 +γ3 ) a 2 u 1 2 1 2 2 a γ2 e2 + u + γ3 e3 + u 2 2
(22)
724
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
ua 2 γ2 ua 2 γ3 2 2 F (v) = a p (γ2 + γ3 ) a u + 1 2 1 2 2 a u γ2 e2 + u + γ3 e3 + u + p 2 2
(23)
a 2 θΓ a 2 (1 − θ) Γ . S (v) = 0 2 a Γ 1H − 2aΨ
(24)
and
Eq. (21) is a system of hyperbolic conservation laws and may be converted to an equivalent set of wave equations that represent nonlinear waves propagating at characteristic speeds, which in turn are functions of the local solutions, and vary in space and time. These characteristic speeds are λ 1 = u − cs ,
λ2 = λ3 = u,
λ4 = u + cs
(25)
where cs is the speed of sound: cs2 = B1 B2 , 1 , B1 = 2 a (γ2 + γ3 )
Rγ3 B2 = a 2 p 1 + C2 γ2 + C3 γ3
! .
(26)
It is also useful to express the equations in terms of the four primitive variables, γ2 , γ3 , u and p. This facilitates the conversion of the boundary conditions to a form that is more convenient for numerical computational purposes. In terms of these variables, the equations are ∂ a 2 γ2 = −d1 + Y1 , ∂t ∂ a 2 γ3 = −d2 + Y2 , ∂t ∂u = −d3 + Y3 , ∂t ∂ a2 p = −d4 + Y4 ∂t
(27)
where a 2 γ2 (Υ1 + Υ4 ) + Υ2 , B2 a 2 γ3 d2 = (Υ1 + Υ4 ) + Υ3 , B s2 B1 d3 = (Υ4 − Υ1 ) , B2 d1 =
d4 = Υ1 + Υ4 ,
(28)
725
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
# ∂ a2 p cs ∂u − , ∂x B1 ∂ x " # ∂ a 2 γ2 a 2 γ2 ∂ a 2 p Υ2 = λ2 − , ∂x B2 ∂x " # ∂ a 2 γ3 a 2 γ3 ∂ a 2 p − , Υ3 = λ3 ∂x B2 ∂x " # λ4 ∂ a 2 p cs ∂ (u) Υ4 = + , 2 ∂x B1 ∂ x λ1 Υ1 = 2
"
(29)
and Y1 = a 2 θΓ , Y2 = a 2 (1 − θ) Γ , uΓ , Y3 = − γ2 + γ3 Y4 =
Rγ3 γ2 C2 + γ3 C3
(30) !
p u2 a 2 Γ p (1 − θ) . a Γ 1H − [θ (C2 − C3 ) + C3 ] − T0 + − 2aΨ + 2 γ3 Rγ3 2
Finally, we denote the temperature within the tube wall as ψ. Inside the tube wall we have the heat diffusion equation in cylindrical coordinates (r, θ, x), 2 ∂ ψ ∂ψ ∂ψ 1 ∂ = Dp r (31) + ∂t r ∂r ∂r ∂x2 where D p = 1.03 × 10−7 m2 /s is the thermal diffusivity of the wall. Let y = r − a0 be the distance measured from the tube wall, where a0 = 1.5 mm denotes the tube radius. Now the length scale that the heat penetrates into the tube is p [y] = D p tr = 2.03 × 10−6 m where tr = 40 µs is the time scale for reaction. We can see that [y] a0 . This allows us to approximate the second term on the right-hand side of (31) by ∂ 2 ψ/∂ y 2 . Also, since the length scale for reaction lr = 0.08 m is much larger than [y], we can neglect the axial conduction term, ∂ 2 ψ/∂ x 2 . The heat diffusion equation (31) thus becomes ∂ψ ∂ 2ψ = Dp 2 , ∂t ∂y
(32)
where ψ = ψ(x, y, t). Note also that the length scale [y] of penetration is a lot less than the thickness h of the tube, since h = 8.5 × 10−4 m. Thus, the far region of the tube wall effectively remains unheated as the detonation wave passes. As a consequence of this, we are able later to reduce the domain of integration for the wall and hence speed up computations. 4. Numerical method 4.1. Initial and boundary conditions Each eigenvalue λi obtained in the previous section represents the characteristic speed at which a particular wave mode propagates. At a point on the boundary, some of the characteristics describe outgoing waves, while some of them describe incoming waves. The behaviour of the outgoing waves is completely determined by data contained
726
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
within and on the boundary of the computational domain. The number of boundary conditions that must be specified at a point on the boundary is equal to the number of incoming waves at that point. We specify the boundary conditions that determine the value of Υi for incoming waves, and compute, from the definition (29), the value of Υi for outgoing waves. At the open end, the gas is free to flow out of the tube. The number and type of boundary conditions required may vary from time to time at this boundary. There are two possible cases in this study. (1) Supersonic outflow: supersonic outflow at the open end is characterised by u > cs , so that all λi are greater than zero. Consequently, we can specify no boundary conditions at all, and the evolution of the flow at the boundary is determined completely by interior data. (2) Subsonic outflow: subsonic outflow is characterised by 0 > u > −cs , for which λ1−3 < 0 and λ4 > 0. Thus, we compute Υ1−3 from their definitions in (29) and specify Υ4 according to some boundary condition. In this study we specify the pressure to be one atmosphere. We achieve this by requiring that ∂(a 2 p)/∂t = 0. This condition is satisfied if Υ4 + Υ1 − Y4 = 0, giving Υ1 = Y4 − Υ4 . The other end of the tube is closed. The initial state has zero speed at the closed √ end, and the boundary conditions are chosen to keep the speed zero at all time. This condition is satisfied if − B1 /B2 (Υ4 − Υ1 ) + Y3 = 0, giving √ Υ1 = B2 /B1 Y3 + Υ4 . We choose the initial conditions to simulate hot combustion gases over a short distance L 1 near the free boundary, i.e. 0 ≤ x ≤ L 1 L. The gas pressure is taken to be constant throughout the tube at the external ambient value. A linear variation of gas temperature is assumed to be present near the open tube end: x , 0 ≤ x ≤ L 1, T0 + (Ti − T0 ) 1 − (33) T (x, 0) = L1 T0 , L1 < x ≤ L , where T0 is the reference temperature (nominally 293 K) and L 1 equals 0.1 m. This suffices to initiate the combustion process inside the tube, provided that the temperature Ti resulting from the ignition source exceeds the ignition criterion for the HMX. For the tube wall, we assume that there is no heat loss to the surrounding medium, so we impose an insulation condition at infinity. At the inner tube wall boundary, we need to specify two boundary conditions. These are continuity of the temperature field ψ (x, 0, t) = T (x, t) , and continuity of heat flux ∂ψ −Ψ = K p ∂ y y=0
(34)
(35)
where K p = D p C p ρ p is the wall conductivity, C p is the specific heat of the tube, and ρ p is density of the tube. As usual, the wall is assumed to be initially at ambient temperature throughout. This may be written as ψ (x, y, 0) = T0 .
(36)
4.2. Numerical solution procedure We use the time splitting method to obtain the numerical solution of Eq. (21). The first step is to remove the inhomogeneous term S(v) for Eq. (21). Thus the system reduces to ∂v ∂F (v) + = 0. ∂t ∂x
(37)
Eq. (37) can be solved using the TVD MacCormack method [16]. Let v¯ be the numerical solution of Eq. (37). Then the system of ordinary differential equations given by dv = S (¯v) , dt
(38)
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
727
together with the reaction Eq. (18), are solved in turn. Eq. (38) is solved using a Runge–Kutta Scheme by Cash and Krap [17], by utilizing the solution v¯ of Eq. (37) to determine the inhomogeneous term S. The heat conduction equation in the wall is parabolic in the time variable. We solve it by sequential application of the method of lines, integrating in any one step from time tn to tn+1 . We define 1y to be the mesh spacing within the wall. Let m 1 denote the number of x-mesh spacings and m 2 denote the number of y-mesh spacings. Then m 1 = L/1x and m 2 ≤ h/1y. We introduce mesh functions defined as the solution of the finite difference equations below and which are related to the exact solution by ψi, j ' ψ(i1x, j1y, t), where 0 ≤ i ≤ m 1 , 0 ≤ j ≤ m 2 . Thus, using standard three-point central differencing, we obtain the following ordinary differential equation (ODE): ψi, j+1 − 2ψi, j + ψi, j−1 dψi, j = . dt (1y)2
(39)
At each time step, we set ψi,0 = T (i1x, tn ) and ψi,m 2 = T0 . Since the heat is only absorbed in a very thin region near the tube wall, we choose m 2 to equal one tenth of h/1y, i.e. we only compute the temperature over one tenth of the tube thickness in our computation. This saves a lot of computational time. Eq. (39) consists of a system of ODEs that are solved using the same Runge–Kutta method. Thus, the inhomogeneous term and the heat loss term are solved together. 5. Numerical results To demonstrate the importance of the tube wall interactions, three sets of results will be presented here. The first set is for results without any tube wall interactions, the second set of results is with wall heat loss only, and the third set is for the complete model. All the solutions are depicted at equal time intervals of 0.1 ms. The tube length is L = 2 m and 801 nodes are used in the calculation. 5.1. Without any tube wall interactions In this section, the numerical results without any tube wall interactions are presented, i.e. the radius a is a constant and Ψ = 0. The graphs of γ1 , γ2 , γ3 , p, T and Mach number as a function of x are shown in Fig. 1. The very last profile shows conditions after the shock has hit the closed end of the tube and reflected. The HMX concentration Fig. 1(a) shows that the reaction zone extends about 0.11 m. The aluminium flake and gas concentration are shown in Fig. 1(b) and (c), respectively. We can see that they have similar features. This is because of our assumption that the gas and aluminium flake move with the same speed. The peak pressure in the pressure Fig. 1(d) is about 40 MPa and the pressure at the open end stays nearly constant over the calculation. Fig. 1(e) depicts the solution for the temperature. The peak temperature is about 5300 K. The Mach number Fig. 1(f) shows that the gas exits at sonic speed and the Mach number near the front is about 0.7. The detonation speed is about 2800 m/s. From our numerical results, the peak pressure is about 40 MPa and the detonation speed is about 2800 m/s. However, if we compare these with information from ICI, we find that our peak pressure is about four times more than the ICI observations. Also, our numerical detonation speed is faster than their experimental observation of 2000 m/s. 5.2. The effect of heat loss to the tube wall In this subsection, we have included the heat loss term in our model. The numerical solutions are shown in Fig. 2. The pressure result in Fig. 2(a) shows a sharp discontinuity. The peak pressure is about 18 MPa, which is a reduction of over 50% in magnitude compared to the result without heat loss. The detonation speed is about 2000 m/s, which is now in good agreement with the ICI experiment results. The gas temperature Fig. 2(b) shows that the peak temperature is about 3200 K, which is about half of the maximum temperature rise without the heat loss effect. The aluminium flake and gas concentration are shown if Fig. 2(c) and (d), respectively. We can see that they have smaller peak values in comparison with the results without heat loss. Fig. 3 shows the temperature inside the wall. As we expected, the heat is only absorbed within a very thin region near the interior tube wall.
728
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
(a) HMX concentration.
(b) Al/Al2 O3 concentration.
(c) Gas concentration.
(d) Pressure.
(e) Gas temperature.
(f) Mach number. Fig. 1. Numerical results without any tube wall interactions at equal time intervals of 0.1 ms.
5.3. Full model results The full model results are shown here. The graph of pressure is shown in Fig. 4(a). We see that the peak pressure is about 13 MPa. This is about a one third reduction compared with the previous result. Physically, this is because the gas now occupies a bigger volume. The corresponding temperature field is shown in Fig. 4(b). There is not much difference between this and the temperature graph without tube expansion from the previous section. This is because we have no energy loss from the tube expansion. The graph of the inner tube radius is shown in Fig. 4(c). We see that it shows a similar profile to the pressure graph. This is because they are connected by the linear relation equation (5). There is about a 14% increase in the tube radius at the peak pressure.
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
(a) Pressure.
(c) Al/Al2 O3 concentration.
729
(b) Gas temperature.
(d) Gas concentration. Fig. 2. Numerical results for heat loss model at equal time intervals of 0.1 ms.
Fig. 3. Tube wall temperature profile at 0.8 ms.
The solutions for the gas concentration and aluminium particles concentration are given in Fig. 4(d) and (e), respectively. We see that, as before, there is a constant state in between the detonation and the rarefaction. Also,
730
(a) Pressure.
(c) Inner tube radius.
(e) Al/Al2 O3 concentration.
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
(b) Gas temperature.
(d) Gas concentration.
(f) Gas speed. Fig. 4. Numerical results for full model at equal time intervals of 0.1 ms.
it is more noticeable than before. The gas speed is shown in Fig. 4(f). The typical gas speed is about 800 m/s after the shock wave. 6. Discussions Our model has a ZND wave structure, that is, a shock followed by a deflagration wave in the HMX, including the tube wall heat loss effect and the tube expansion effect. The numerical results of our model compare favourably with experiments performed by ICI. However, there are some aspects of the results that we need to discuss further. If the pressure results in Fig. 4(a) are examined carefully, a pulsation detonation pressure profile is observed. Fig. 5 shows the maximum detonation pressure against time, which clearly reveals that the model is unstable. The frequency
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
731
Fig. 5. Maximum detonation pressure against time (dashed line, 801-nodes; solid line 1601-nodes).
of the oscillation remains approximately constant when we enlarge the number of nodes. Thus we believe that the instability is not caused by the methodology. As we mentioned at the beginning of this paper, the ZND wave structure is usually unstable [13] and our numerical calculations clearly show a one-dimensional instability. Another interesting feature is the constant state in between the reaction zone and the head of the rarefaction wave: in the first set of results the double heads emerge in the solution. This is due to the final state of the detonation moving too fast such that it runs away from the following rarefaction wave. We can see from the pressure graph Fig. 1(d) that the first head begins to develop at 0.4 ms and its strength begins to increase, whereas the second head strength starts to decrease. Also, the distance between the two heads increases and the detonation moves further away from the rarefaction wave. The double-head feature disappeared in the second set of results, but there is a distinguishable region in between the detonation and the head of the rarefaction wave. The wall heat loss effect caused some energy in the gas to be transferred to the tube wall. This reduces the strength of the detonation front and the rarefaction wave. However, heat loss has much effect on the rarefaction wave and reduces the wave propagation speed relative to the detonation. Hence the second head vanishes in the solution and we can see a constant state. This constant state becomes more noticeable when we add the tube expansion effect. The speed of propagation of the rarefaction waves in between the second and third set of results is about the same. But this time the tube expansion significantly reduces the detonation head strength, such that we see the constant state more clearly, although there are some wrinkles along the constant state. We believe that this is caused by the instability of the model. This constant state tells us that our detonation speed is faster than the CJ speed and our final result has a weak solution. 7. Conclusions In this paper we have formulated a one-dimensional model in which the detonation wave is propagated in detonation transmission tubing and where the reaction rate is described by a single-step kinetic equation. We then went on to describe the characteristic and the numerical method of the model equation. The heat loss effect and the tube expansion effect have also been considered. The final results obtained are in good agreement with the experimental data. The results show that the heat loss and the tube expansion have a major effect on the structure of the detonation wave. The long reaction zone results in a non-adiabatic combustion system. The elastic material of tubing expands under high interior pressure. Also, the constant state in between the detonation and the rarefaction wave tell us that we have a weak detonation. An instability feature causes an unstable (steady oscillation) detonation. Acknowledgements The author would like to express his gratitude to ICI for their generous financial support and to Dr A. Janes for his valuable help during this project. The author would also like to acknowledge the highly useful discussions about this problem with Dr. D. Sutton and Dr. P. Lynch of ICI Explosives. Appendix In this appendix, we derive and discuss the various physical scalings that are quoted in the text and used to justify the various assumptions made.
732
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
One physical quantity needed is the time scale of reaction. We were provided with the results of some experiments on DTT performed by ICI Explosives Division which enabled us to make an estimate of this. In one of these experiments, the light emitted by the detonation was measured as it passed a fixed point of the tube. It seems natural to assume that this light is generated by the reaction zone, and that the light zone τ is the required reaction time. The reaction constant φ in Eq. (18) can then be chosen to generate such a reaction by choosing φ = 3(γ10 )1/3 /τ , where γ10 = 11.09 kg/m3 is the initial concentration of HMX. The observed light period is about τ = 100 µs. The use of this value in our numerical calculations, however, predicted a detonation speed significantly below the experimentally observed speed of 2 km/s. Also, the intensity of the light emission profile does not fit the theoretical profile well in the early stages. There is a 15 µs lag before the maximum light output was achieved, whereas the theoretical profile has its maximum attained immediately. Such disagreements are not entirely unexpected, given the highly simplified assumptions made about the kinetics and the geometry. Eventually, we elected to choose a value for τ that gave good agreement with the observed detonation speed, rather than with the light emission, and the value satisfying this criterion, τ = 40 µs, was found by numerical experiments. One of the assumptions made was that the speed of a typical aluminium particle was the same as that of the gas flow. A typical Reynolds number for an aluminium flake is ρulal = 1691 µ
Ral =
where u ∼ 800 m/s is a typical gas speed (taken from the final numerical results), ρ ∼ 16 kg/m3 is a typical gas density after detonation (also from the numerical results), µ ∼ 7.57 × 10−5 kg/ms is the dynamic gas viscosity (of air) at 2500 K, and lal = 10 µm is the relevant length scale for gas flow past a thin plate. Since the Reynolds number Ral is large, the motion of the aluminium obeys inviscid dynamics. In inviscid flows the acceleration of the particle is caused by the pressure difference on the forward and backward faces of the particle. This pressure difference has the order of magnitude 21 ρu 2 , so let us take a mean pressure difference as 1 1 2 1p = ρu . 2 2 For an aluminium flake, this acts over an area of lal2 /2 (because of the aspect ratio effect). Therefore, the typical force is Fal =
lal2 ρu 2 . 8
By Newton’s second law, Mal
du al = Fal , dt
where Mal is the mass of each aluminium particle. So the time 1tal needed to achieve a speed change 1u al = u is 1tal =
Mal u 8ρal dal = = 1.68 × 10−7 s, Fal ρu
where dal = 0.1 µm is the thickness of the aluminium flake and ρal is the density of aluminium. We can see that the time scale is fast compared with the reaction time scale tr = 40 µs. Hence we assume that aluminium particles move with the speed of the gas. A similar calculation for HMX, which has density ρHMX = 1900 kg/m3 , gives 1tHMX = 2.71 × 10−5 s. We have assumed that the HMX is stationary, but clearly this not as good an approximation as the last one. The particles would be accelerated to the gas speed over a time comparable to the reaction time, and ideally the model should be modified. However, the chemical reaction is proportional to the surface area, so most of the reaction will take place during the early stages of the acceleration process when the speed that has been achieved is small. Hence we expect the effect of the motion to be relatively minor, and so neglect it.
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
733
The time scale to heat up a HMX particle by diffusion is 2 1 2 lHMX = 4.0 × 10−3 s, [t] = DHMX where DHMX = 1.0 × 10−7 m2 /s is the diffusivity of the HMX. This is slow compared to the time scale for the reaction. We can conclude from this that the bulk of the particle does not have time to heat up by diffusion, and hence that it does not react as a uniform whole, and must therefore react by conversion of the surface of HMX to gas. Heat diffusion takes place, but only in a thin layer next to the newly exposed surface. This is rapid because of the large temperature gradient. Diffusion into the bulk is negligible. We thus assume in the model that the rate of reaction Γ is proportional to the surface area of HMX remaining, and that the temperature of the unburnt HMX remains constant during the reaction. Conversely, the time scale taken to heat up the aluminium flake of thickness dal is 2 1 2 lal = 2.56 × 10−11 s, [t] = Dal where Dal = 9.76 × 10−5 m2 /s is the diffusivity of the aluminium. Now this value is very short compared with our reaction time scale, so we assume that the temperature of the aluminium flake is the same as that of the gas. From thermodynamic principles, the oxygen reacts with the metal preferentially to form metal oxide, then any remaining oxygen reacts with carbon to form carbon monoxide, and then the remnant oxygen reacts with hydrogen to form water. According to this rule, we can write the overall chemical formula per unit mole of HMX as Eq. (1). We can see from Eq. (1) that the reaction product is dominated by carbon monoxide and molecular nitrogen. We treat the gases as a single gas product, taking the molecular weight of the hypothetical gas to be the mole average of the gas products. The average molecular weight W = 22.99 g/mol is used here. Since the aluminium flake is heated instantaneously, we assume that the oxygen initially present in the air reacts immediately with aluminium to form Al2 O3 , such that no free oxygen remains in the air after the passage of the shock. The chemical formula for this reaction is 0.954Al + 0.155O2 → 0.103Al2 O3 + 0.747Al. We see that only about 22% of that aluminium becomes Al2 O3 from this reaction. To avoid working with a two-reaction kinetic scheme, we decided to omit this reaction. For all the aluminium to become Al2 O3 , we need a further 0.56 mole of oxygen from the 1 mole of HMX. The chemical formula is 0.747Al + 0.560O2 → 0.375Al2 O3 . weight of oxygen So the mass ratio θ = 0.56×molecular = 0.061. We can see that only a very small fraction of gas is molecular weight of HMX absorbed by aluminium. Most of the HMX becomes gas.
References [1] M. Berthelot, P. Vielle, On the velocity of propagation of explosive processes in gases, C. R. Hebd. Sceances Acad. Sci. 93 (1881) 18–21. [2] M. Berthelot, P. Vielle, On explosive waves, C. R. Hebd. Sceances Acad. Sci. 94 (1882) 882. [3] E. Mallard, H. Le Chatelier, On the propagation velocity of burning in gaseous explosive mixtures, C. R. Hebd. Sceances Acad. Sci. 93 (1881) 145–148. [4] D.L. Chapman, On the rate of explosive in gases, Philos. Mag. 47 (1899) 90–104. [5] E. Jouguet, On the propagation of chemical reactions in gases, J. de Math. Pures et Appl. 1 (1905) 347–425. Continued in 2 (1906) 5 – 85. [6] E. Jouguet, Mecanique des Explosifs, Octave Doin et Fils, Paris, 1917. [7] Y.B. Zeldovich, On the theory of the propagation of detonation in gaseous system, Zh. Eksp. Teor. Fiz. 10 (1940) 542–568. English translation: NACA TM1261, 1960. [8] J. von Neumann, Theory of detonation wave, Technical Report, OSRD Report No. 549, Paris, 1917. [9] W. Doering, On detonation processes in gases, Ann. Phys. 43 (1943) 421–436. [10] A.K. Oppenheim, J. Rosciszewski, Determination of the detonation wave structure, in: Proc. 9th Symp. Int. on Combustion, 1963, pp. 424–441.
734
D.K.L. Tsang / Mathematical and Computer Modelling 44 (2006) 717–734
[11] D.H. Edwards, G.T. Williams, J.C. Breeze, Pressure and velocity measurements on detonation waves in hydrogen–oxygen mixtures, J. Fluid Mech. 6 (1959) 497–517. [12] J.H.S. Lee, Dynamic parameters of gaseous detonation, Ann. Rev. Fluid Mech. 16 (1984) 311–336. [13] G.E. Abouseif, T.Y. Toong, Theory of unstable one-dimensional detonations, Combust. Flame 45 (1982) 67–94. [14] F. Zhang, H. Gronig, Spin detonation in reactive particles-oxidizing gas flow, Phys. Fluids A 3 (1991) 1983–1990. [15] T.B. Brill, H. Arisawa, P.J. Brush, P.E. Gongwer, G.K. Williams, Surface chemistry of burning explosives and propellants, J. Phys. Chem. 99 (1995) 1384–1392. [16] S.F. Davis, TVD finite difference schemes and artificial viscosity, Technical Report, NASA CR 172373, 1984. [17] J.R. Cash, A.H. Karp, A variable order Runge–Kutta Method for initial value problems with rapidly varying right-hand sides, ACM Trans. Math. Software, 16 201–222.