A stability analysis is carried out taking into account slightly porous walls in an idealized detonation confined to a circular pipe. The analysis is carried out using the normal-mode approach and corrections are obtained to the underlying impenetrable wall case results to account for the effect of the slight porosity. The porous walls are modelled by an acoustic boundary condition for the perturbations linking the normal velocity and the pressure components and thus replacing the conventional no-penetration boundary condition at the wall. This new boundary condition necessarily complicates the derivation of the stability problem with respect to the impenetrable wall case. However, exploiting the expressly slight porosity, the modified temporal stability can be determined as a two-point boundary value problem similar to the case of a non-porous wall.
The normal-mode analysis was the pre-eminent feature of conventional hydrodynamic stability theory for a number of years before it was recognized that the receptivity problem is also an important element of the flow response [1,2]. The analysis requires the solution of partial differential equations (PDEs), and a possible formulation is the initial-value problem (IVP) where one specifies the initial profile of the perturbation flow. Using the Laplace transform with respect to time leads to a simplified set of equations and a formal solution in the transformed space. A detailed analysis of the inverse Laplace transform can provide expansion of the solution into normal modes of discrete and continuous spectra. The weights of these discrete modes represent the receptivity of the flow to the particular initial perturbation field. This type of analysis in detonation stability was not completed in Erpenbeck's original work [3,4] since the main interest at that time was in the singularities associated with the unstable modes. In the unbounded case, Tumin [5; 6] completed the relevant receptivity analysis for one-dimensional perturbations of the idealized detonation model, and Chiquete et al. [7,8] extended this analysis for three-dimensional perturbations of the same flow for overdriven and explicitly Chapman–Jouguet (CJ) detonation.
Experimental studies of detonations often take place within the confining pipe geometry. Experimental studies involving detonations in porous-walled tubes have been carried out as well, notably by Teodorczyk & Lee . They found that the porous wall is crucial in the attenuation of transverse waves, which have been observed to contribute to the persistence near the detonability limits of the propagating detonation wave.
Radulescu & Lee  provided further insight: detonations in highly unstable mixtures confined by porous walls fail primarily through the corresponding attenuation of these transverse waves. Conversely, detonations in weakly unstable mixtures fail owing to global curvature effects caused by the mass divergence into the porous wall. The attenuation of the transverse waves is observed in numerical simulation as well . The primary motivation in understanding the attenuating properties of these porous-walled pipes is their observed ability to quench unwanted detonations that can develop in explosive gas mixtures that occur in certain chemical industrial settings .
To our knowledge, there has been no stability analysis carried out within the porous-walled tube configuration. The effect of a porous wall on the stability problem for detonations will therefore be explored in the present work within the linearized stability framework. The porosity is modelled through use of the so-called acoustic boundary condition where the normal velocity component of the perturbation flow is no longer zero at the boundary but proportional to the pressure component. As a first-order approximation, the analysis is restricted to the case where the pipe walls are only slightly porous in order to maintain the basic mathematical structure of the limiting case of the non-porous wall and facilitate a comparison between the two cases. Crucially, the approximation allows the definition of the stability problem as a two-point boundary value problem as in the impenetrable wall case and thus enables the use of the same numerical methods with only slight modification.
2. Governing equations
The full set of governing equations (in the laboratory frame) for idealized one-reaction detonation are the Euler equations supplemented by an equation governing the reaction kinetics, 2.1a 2.1b 2.1c 2.1d where D/Dt=∂/∂t+ul⋅∇l and superscript l denotes the laboratory frame, is the reaction heat release and r is the rate of the reaction. Here, u=(u,w,v)T, i.e. radial, azimuthal and axial velocity. The usual symbols denote pressure (p), specific density (ρ), reaction progress variable (λ) and specific internal energy (e). The equations and subsequent dimensional scales are introduced following Kasimov & Stewart . For an ideal poly-tropic gas, 2.2 where λ represents the reaction progress variable, with λ=1 denoting the unburned state and, naturally, λ=0 denoting a totally burned state. The ideal equation of state is used here, which leads to the sound speed c, density, specific heat ratio γ and pressure having the relation c2=γp/ρ and the Arrhenius reaction rate is used as well, 2.3 where is the temperature, is an activation temperature, μ is the mean molecular weight of the gas and finally is the specific gas constant. The Rankine–Hugoniot jump conditions hold at the shock front surface.
One can then posit that there is a steady one-dimensional in the frame of reference moving with the shock. This flow is commonly referred to as the Zel'dovichvon NeumannDoering (ZND) idealized model (see Fickett & Davis  or Williams ). Therefore, the flow supposes that there is no steady-state velocity in the radial or azimuthal directions. Additionally, the superscript ‘*’ is used to denote this base flow.
Note that the scales are also consistent with Kasimov & Stewart , i.e. half reaction length and post-shock values of density, pressure and sound speed. There are four dimensionless parameters that characterize the ZND flow summarized here:
— : overdrive factor measuring the extent by which the detonation speed exceeds the minimally allowed Chapman–Jouguet velocity .
— : activation energy scaled with the quiescent speed of sound squared.
— : heat release scaled with the quiescent gas sound speed squared as in the above.
— γ: specific heat ratio.
(a) Stability analysis
In the current approach, the normal-mode ansatz is used. A detonation propagates to the right in one direction (zl), and its position is modified slightly by ψ′, which is a small perturbation. One proposes a change of coordinate to the moving perturbed shock, 2.4 where Ds is the detonation velocity scaled with the post-shock sound speed. One can also choose to measure particle velocity with respect to the unperturbed shock, 2.5 Additionally, one supposes that the ZND base flow is modified by small perturbations which fully depend on the axial coordinate z, radial coordinate r, azimuthal angle ϕ and time t. For example, one holds that the density is expanded as 2.6 The full ZND flow is given in vector form as 2.7 and the solution to the full Euler equations (2.1) is expanded via the standard linearizing procedure, i.e. 2.8 The perturbation vector q is considered to be very small and, as a result of ignoring the higher order terms, one obtains the system of PDEs, 2.9 where q′=(ρ′,u′,w′,v′,p′,λ′)T is the relevant vector of perturbations, and the coefficient matrices and vectors are functions of the base flow. The appropriate boundary condition on the shock front (z=0+) is derived from the linearized Rankine–Hugoniot jump conditions, 2.10 At the end of the reaction zone, a spatial boundedness constraint is imposed, 2.11 This is formally required as but, practically, only a few reaction lengths are sufficient to reach the end of the reaction zone. The normal-mode approach to determine stability relies on the following ansatz: 2.12 As a result of this choice, one encapsulates the temporal instability of the flow in the real part of τ and its frequency in the imaginary part of τ as well. Under this assumption, the task becomes to obtain a reduced eigenvalue problem that will determine τ for the particular detonation in question. Ideally, one would hope for the governing equations of this eigenvalue problem to be a set of ordinary differential equations (ODEs). To that end, the procedure continues via a Fourier series expansion in the azimuthal variable in order to eliminate the explicit dependence on ϕ, 2.13 These transformations lead to a system of PDEs in r and z but now lacking derivatives in t and ϕ for the reduced eigenfunctions zn and ψn, 2.14 The transformed boundary conditions in the z direction are, specifically, 2.15
(b) The eigenfunctions in the non-porous wall case
The problem formulation is incomplete until a boundary condition on the wall is defined. The governing equations depend on both r and z and, formally, the resulting eigenvalue problem is paired to two-dimensional eigenfunctions. Ideally, a further ‘separation of variable’ scheme could produce a reduced set of ODEs to define the eigenvalue problem rather than the more complex set of PDEs in (2.14). This separation scheme should of course depend on the boundary condition in r. Specifically for the conventional case, there is a no-penetration boundary condition in the radial velocity on the walls (r=a), 2.16 The boundary condition above specifically allows the definition of the following separation scheme based on Bessel functions and its derivatives (where Jn refers to the Bessel function of the first kind of order n; table 1)Q5 found by Kasimov & Stewart , i.e. 2.17 and 2.18 The eigenvalues knj come about from the boundary condition at r=a, 2.19 In our previous work in the linear stability of the ZND flow in the unbounded case, the wavenumber k in the transverse direction belonged to a real, continuous set. Conversely, in the restricted finite domain of a pipe geometry, the parameter becomes a discrete number and is obtained from the no-penetration wall condition. However, note that Shalaev & Tumin  show that the stability of the confined case is equivalent to the unbounded case via a transformation of variables. Additionally, they find that the unstable modes carry no axial vorticity by definition, pointing to the importance of the continuous spectrum that arises in the IVP approach but not within the normal-mode analysis of the discrete spectrum.
The set of radial wavenumbers knj introduces the orthogonality relation, 2.20 and this leads to a system of ODEs in z for the quantities . Each of these can then be placed into a vector pj and one can then obtain the two-point boundary value problem, 2.21 where T and F depend explicitly on the base flow quantities, knj and n. These eigenfunctions pj and their corresponding eigenvalues τj determine the stability of the ZND base flow subject to their compliance with the previously expressed boundary conditions. This particular problem is solved in Kasimov & Stewart  and is also treated via the IVP approach in Shalaev & Tumin .
3. The porous wall and the acoustic boundary condition
The focus of the present work is to model the effect of a porous wall on detonations within the framework of the linearized stability analysis. The model involves incorporating an acoustic absorbing boundary condition on the wall (i.e. r=a) in lieu of the no-penetration condition used up to this point in order to model surfaces that are neither rigid nor impenetrable. The boundary condition is built on the assumption that the perturbation pressure should track linearly with the normal velocity at the surface in question . The specific acoustic impedance is the quantity that links these two components, and its introduction to acoustics dates to Webster  in 1914. Therefore, each component perturbation of a prescribed frequency is uncoupled from the rest via this linear relationship , i.e. one has for each of these temporal modes the following: 3.1
where and denote the pressure and normal velocity component mode for each defined time frequency ω at the normal surface, S0 is the normal surface in question and Zs(ω) is the frequency-dependent specific impedance (complex number in general). The formulation of the acoustically absorbing boundary condition similar to (3.1) first appeared in 1921 and it is due to Kennelly & Kurokawa .
Matching the acoustic boundary condition that results from the analysis of the porous layer to the outer inviscid flow where the ZND flow applies remains to be examined. In problems involving acoustic perturbation propagation in ducts, the conventional approach has been to model the boundary layer as a thin vortex sheet and where appropriate jump conditions are found across the thin interface. This condition was originally used and derived by Taylor  in 1979 and, later, more formally derived by Myers . It involves matching the perturbation normal velocity to the normal velocity of the surface at each point and thereby ensuring continuity of particle displacement. One obtains then the boundary condition of the form (following the notation used by Farassat & Dunn ), 3.2 where u1 is the normal velocity perturbation in the inviscid flow and u0 is the base flow in the same region. The surface S(t) is assumed to vary weakly away from a mean surface S0, where g(x,t) is the perturbation to the normal displacement of the surface and for the current case and n0 is the normal vector to the mean surface. In the current work, the normal direction to the surface is parallel to the radial coordinate and the mean surface lies at r=a. In the one-dimensional ZND flow, there is only a base velocity in the axial direction, 3.3 and it is held that the normal velocity ‘below’ the thin boundary layer surface tracks with the perturbation pressure according to (3.1), i.e. in the porous layer, . The time dependence of each of the relevant quantities is assumed proportional to . Since , it follows that . Continuity of pressure across the interface is also held so that the perturbation pressure in the porous and inviscid flow are equal. In the shock frame of reference then, 3.4 The impedance Z will be slowly varying in the axial variable z and therefore the term stemming from its gradient is negligible for the current case. The boundary condition for the inviscid perturbations is therefore, 3.5 The specific geometry of the walls will affect the impedance Z and, additionally, the eigenvalue (i.e. τ) of the normal-mode perturbation. Finally, note that in the case of infinite impedance the impenetrable wall case boundary condition is recovered.
This new acoustic boundary condition unfortunately complicates the mathematics considerably. As a result, the stability is formally determined by solution of a PDE system in r,z for two-dimensional eigenfunctions and eigenvalue τj, all parametrized by n. Clearly, both the computational complexity and cost of the numerical implementation for the stability problem solution increase with respect to the impenetrable wall case. However, as a first step, an approximate analysis of the stability problem was undertaken in order to ascertain the effect of slight porosity on the stability problem.
In addition to the porosity effect at the boundary wall, the slight relaxation of the no-penetration condition can introduce a further effect upon the governing equations for the perturbations as well. The small mass divergence into the wall precipitated by the acoustic boundary condition leads to slight curvature of the leading shock front . Nevertheless, this leads to the modification of the linearized governing equations for the interior flow perturbations at , where κ is inversely proportional to the radius of curvature. In the work of Radulescu & Lee , the mass divergence introduced by the porous walls is related to the corresponding induced curvature measure κ as 3.6 where β is the open area ratio of the porous wall, which in this study is considered to be small and inversely proportional to the acoustic impedance, and d is the measure of the height of the channel. This estimate is derived based on the assumption of a weakly curved ZND shock front profile and relating the radial gradient of the normal velocity component to the relative rate of area increase in the reaction zone . This relation reveals that curvature effects upon the linear stability analysis in the porous wall case are slight when the height of the detonation channel is large compared with the reaction zone length. In the present work, the effects of curvature are therefore neglected and so the corresponding results should be interpreted as the first approximation in the large though finite channel height limit.
(a) Mathematical approach to the stability question
The overriding philosophy to determine the stability of the reduced problem for a slightly porous wall is through an eigenfunction expansion in r (as described in general, for example, by Zauderer ). In the current case then, the jth normal mode is expanded in the radial variable using the eigenfunctions obtained from the non-porous case, 3.7 and 3.8 In order to obtain the weight, , , etc., for each radial eigenfunction, one incorporates a type of integral inner product over the radial domain. This inner product is incorporated equation by equation by integrating against the appropriate radial eigenfunction. For example, one has for the continuity equation 3.9 parametrized by n, the azimuthal wavenumber, and obtained subsequent to the implementation of the exponential time dependence prescribed by the normal-mode approach. One proceeds by multiplying through by the function Jn(knjr)r and integrating from r∈[0,a], 3.10 Integration by parts explicitly introduces the boundary condition at r=a, 3.11 Applying the acoustic boundary condition, 3.12 Clearly, as , the equivalent equation to the solid-wall case is recovered and the eigenfunction expansion collapses to a single element of the radial eigenfunction set; in this case, the jth such element.
The weight functions such as are functions only of z. Note that the eigenfunction expansion leads to an infinite sum over the radial eigenfunction indexed by m because replacing in equation (3.12) with the eigenfunction expansion leads to the following: 3.13 The sum is well defined and converges uniformly and absolutely as long as is smooth (minimally: the first two derivatives in r must exist according to Agrawal & Patel ). The integrals in equation (3.12) disappear owing to the orthogonality of the underlying radial eigenfunctions and so one obtains equations that depend only on z.
The derivations of the other relevant conservation equations proceed similarly and appear completely in the study of Chiquete . The salient aspect is that the evolution of the jth unstable normal mode will depend directly on all the radial eigenfunctions from the non-porous case. For the case where 1/Z is small, this can be improved upon as follows. As , one expects that the eigenfunction expansion reverts to the simple form in the non-porous case. Therefore, 3.14 This in effect forces the eigenfunction expansion weights different from j to be or, more generally, vanish in the limit. The input from the other radial eigenfunctions is thus a higher order effect within the governing equations, 3.15 Dropping the higher order terms leads to a system of ODEs that can be written as 3.16 where the eigenfunction vector . The boundary conditions are correspondingly obtained as a matrix-vector system, 3.17 The stability problem is thus a one-dimensional two-point boundary value problem in z as in the non-porous wall case. The approximate approach considered here is illustrated in the study of Chiquete  via a model problem involving the wave equation on a finite domain.
(b) Model for the acoustic impedance
In the current work, a relatively thin porous layer lines the walls of the cylindrical tube. The schematic of the porous wall layer is from a top view (figure 1). For the specific case presented here, one imagines a layer containing a series of cylindrical cavities of a specified radius separated by a certain distance . Note that other cavity cross sections could be used as well. With the geometry of the porous layer defined, the remaining task is to determine the nature of the specific impedance Z produced by this specific layer, which then determines the modified boundary condition, 3.18 The specific impedance Z depends on the geometry, characteristics of the ambient base flow and, crucially, on the time dependence encapsulated in τj of the perturbation itself. Stinson & Champoux  obtained this quantity for the configuration above for each specific cavity and also as a bulk material. Their method relies on solving the continuum fluid equations in each cavity for perturbations to a steady state. In the case that the cavities are indeed of the order of the molecule mean-free path, Kozlov et al.  provide the same quantities for the non-zero Knudsen number. In the present work, the size of the cavities is assumed to be within the continuum flow regime.
The method of Stinson & Champoux  involves solving perturbation equations in the cavity assuming a perfect gas with a constant equilibrium flow and with finite dynamic viscosity. The solution for the perturbations in this restricted domain is used to obtain the constant of proportionality between the normal velocity and pressure perturbations at the surface. The result for each cavity was then integrated in order to obtain the effective specific acoustic impedance for the surface as a whole.
In the specific case considered here, the radius of the pipe is considered to be much larger than the characteristic size of the porous-walled holes yet large enough that the continuum equations hold in these cavities. Using the appropriate characteristic scaling from the ZND detonation flow (i.e. quantities directly behind the detonation flow among others) and based on the analysis described to this point, one obtains 3.19 with 3.20 The impedance is illustrated in figure 2. Though the equilibrium base flow varies through the reaction zone, one assumes a constant equilibrium state in each cavity assuming that the variation of the base flow is large compared with the size of each cavity. In terms of the quantities that appear in the formulae, the J0 and J2 refer to the zero- and second-order Bessel functions of the first kind, Pr is the Prandtl number and and are the dimensional quantities representing the radius of the cavity and the separation of the cavities, respectively. Additionally, np is a measure of the porosity or the ratio of the cavity volume to the overall volume of the layer. The constant ξ depends on the dimensional detonation quantities and , and (respectively, density, speed of sound directly behind the shock and the characteristic reaction zone length). The constant is the characteristic size of the holes in the porous media that make up the pipe's walls. Furthermore, is the dynamic viscosity of the gas. Note that the porosity parameter np basically constrains the magnitude of the overall effect of the acoustic boundary condition because and 1/Z∝np. The main goal of the present study was to observe the effect of this boundary condition on the growth rate of the perturbations as the porosity is increased from zero, i.e. the non-porous wall case.
4. Numerical results
(a) The method
The numerical method for obtaining the evolution of the growth rate and frequency of the perturbations comprises several steps which are detailed below. The numerical methods employed to find the eigenvalues are implemented in Fortran using the numerical libraries LAPACK and SLATEC.
— A multi-domain spectral collocation method [27,28] is used to obtain the entire spectrum for the related one-dimensional perturbation system without the porous wall (i.e. np=0). The method is not appropriate for higher dimensional perturbation problems owing to the requirement of a linear boundary condition in τj at the end of the reaction zone. However, owing to the spectral expansion, the entire spectrum is determined at once.
— These seed values for the τj are used as initial guesses for the Newton–Rhapson shooting method used by Lee & Stewart , i.e. one starts from the known values at the shock front and shoots towards the end of the reaction zone using the boundary condition as a ‘dispersion relation’.
— As np is increased, this Newton–Rhapson method is used to track the resulting eigenvalue and eigenfunctions. Note that as this parameter is increased the inverse specific impedance, 1/Z, is increased and thus so is the effect of the porosity.
The numerical parameters were varied to achieve convergence in the following results.
There is broad scope for the choosing of specific cases to test the approach. Explosives can be broadly characterized as ‘highly’ and ‘weakly’ unstable depending on the amount of dilution of the reactive elements. The following examples are obtained with help from data obtained from Austin . Additional results appear in the study of Chiquete . However, the physical values of the constants should still be considered largely approximate.
(i) Weakly unstable case
This particular case is aimed at testing a weakly unstable case, and the parameters were chosen to be consistent with an argon diluted mixture with formula 2H2+O2+12Ar. There is only one unstable mode, τ=0.1587284+0.5068546i, and, for the case at hand, the stability of this three-dimensional unstable mode is tested with k01=0.5 and f=1.05 (slightly overdriven). Additionally, the ZND flow parameters are γ=1.6, Q=24.2 and E=77.25, and the radius of the tube is considered to be 7.663 times the reaction zone length (table 2). The evolution as np is increased from 0 to approximately 3×10−4 is shown in figure 3. The eigenfunctions are shown in figure 4.
Different Bessel function orders n and wavenumbers knm affect the growth rate term in the same qualitative manner as in the n=0 case (figure 5).
(ii) Highly unstable case
The eigenfunctions, frequencies and the growth rates were obtained for parameters consistent with the mixture C2H4+3O2+10.5 N2 for varying degrees of porosity. This case has more unstable modes and is of a higher growth rate (figure 6) than the weakly stable case presented previously. Each mode is labelled in order of increasing frequency and note that mode 2 is the unstable element with the largest growth rate in this case. The porosity is varied between 0 and about 3×10−4 for each case. The growth rates are given for the non-zero frequency modes shown in figures 7 and 8. The corresponding eigenfunction are shown in figures 9–12.
An extension of the circular pipe stability analysis is considered where the wall is allowed to take on a slight porosity. The effect is modelled as an acoustic boundary condition for transverse first-order perturbations to the idealized ZND flow in the pipe. A model for the acoustic impedance is incorporated to determine the impedance which is then, crucially, a function of the perturbation frequency. Modified equations are then obtained incorporating the modified boundary condition, and the perturbation eigenfunctions and temporal eigenvalue containing the growth rate and frequency are sought. The method of solution is a radial expansion using the Bessel function-derived eigenfunctions of the non-porous wall case. The leading order effect of the slight porosity is retained in the governing equations, thus simplifying the problem. A modified two-point boundary value problem is then obtained to determine the corresponding stability information. This allows the use of the numerical method used to determine stability in the more conventional case.
Two different detonable gases are approximated via the ZND one-dimensional model, one representing a highly unstable nitrogen gas mixture and the other a weakly unstable argon gas mixture. It is found that the slight porosity attenuates or dampens the growth rate of perturbations for different radius, and eigenfunction wavenumbers in both the radial and azimuthal directions. Also, the temporal frequency of the perturbations is decreased and therefore reduces the azimuthal velocity of the transverse waves at the wall. The reduction in the growth rate of the transverse perturbations seems consistent with experimental evidence of transverse wave attenuation. However, it should be stressed that the effects modelled here correspond to the initial, linear stages of the transverse perturbation's growth.
The clear extension to this work involves the more general case where the acoustic impedance is not constrained to be large in order to facilitate the definition of the conventional eigenvalue problem. To this end, recent advances  in defining linear instability in terms of PDEs for eigenfunctions of two or more dimensions provide a template for treating the current porous wall configuration. The method involves discretization of the partial derivatives that appear in the governing equations and thereafter the solution of large systems of linear equations primarily through Krylov subspace iterative methods. The approach has proved effective for determining instability for problems where the base flow is two dimensional or where boundary conditions for the perturbations are inhomogeneous in one or more directions , necessarily complicating any attempt in defining the linear instability problem as a two-point one-dimensional boundary value problem. The case where complex boundary conditions are essential is of special relevance to the current work owing to the acoustic boundary condition. Therefore, applying the PDE approach outlined in the study of Theofilis  to stability of detonation in porous walls should be efficient.
The authors thank Prof. Edward Kerschen for useful discussion regarding the acoustic boundary condition. The work produced here was made possible with funding for C.C. from the National Science Foundation Vertical InteGration in Research and Education (VIGRE) programme no. 0602173.
(a) The matrices for the PDE system in (2.9)
A1 A2 and A3 and the vectors: A4
(b) The matrices and vectors pertaining to the ODE system in (3.16)
One contribution of 12 to a Theme Issue ‘The physics, chemistry and dynamics of explosions’.
- This journal is © 2012 The Royal Society