1. Introduction
A combustion wave in which exothermic energy release contributes to driving a leading shock is a detonation. Much of detonation modelling has neglected diffusive processes, as their effects are thought to be small in comparison with advection and reaction (Shepherd Reference Shepherd2009), cf. Fickett & Davis (Reference Fickett and Davis1979), Oran et al. (Reference Oran, Weber, Lefebvre and Anderson1998), Hu et al. (Reference Hu, Zhang, Khoo and Jiang2004, Reference Hu, Zhang, Khoo and Jiang2005), Tsuboi, Eto & Hayashi (Reference Tsuboi, Eto and Hayashi2007), Deiterding (Reference Deiterding2009) and Taylor et al. (Reference Taylor, Kessler, Gamezo and Oran2012). However, Singh, Powers & Paolucci (Reference Singh, Powers and Paolucci1999) and Powers (Reference Powers2006), in a two-dimensional study of detonation patterns using a one-step kinetics model, demonstrated that for the reactive Euler equations detonation patterns were grid-dependent. Moreover, the patterns using the reactive Navier–Stokes equations reduced to a grid-independent dissipative structure. In addition, in regions that require high resolution for a one-step detonation in a channel, diffusion played an important role (Mazaheri, Mahmoudi & Radulescu Reference Mazaheri, Mahmoudi and Radulescu2012). Furthermore, it was demonstrated by Romick, Aslam & Powers (Reference Romick, Aslam and Powers2012) that for the one-step kinetics model the long time behaviour of the detonation is affected by the addition of viscosity; these viscous effects delay the transition to instability, and in the regime where multiple frequency oscillations exist, viscosity can play a dramatic role. Moreover, Al-Khateeb, Powers & Paolucci (Reference Al-Khateeb, Powers and Paolucci2013) suggested that hydrogen–air mixtures have reaction length scales present which have time scales associated with them over which both chemistry and diffusion can be important. This suggests that viscosity has a role to play in the dynamics of detonations; however, the role of viscosity in pulsating detonations of a detailed kinetics model has not been quantified in a detailed manner.
In detonations propagating into mixtures at atmospheric pressure, there are reaction dynamics and gradients in flow variables evolving on scales near a micrometre or smaller (Powers & Paolucci Reference Powers and Paolucci2005). For a steady, travelling, inviscid, Chapman–Jouguet (CJ) hydrogen–air detonation propagating into ambient atmospheric conditions, the reaction length scales range from
$10^{-5}$
to
$10^{0}~\text{cm}$
at equilibrium, and near the detonation front this breadth of scales can be even wider. These length scales are a manifestation of the averaged continuum representation of the collisions that occur between molecules, which have an underlying principal length scale of the mean free path (Al-Khateeb, Powers & Paolucci Reference Al-Khateeb, Powers and Paolucci2010). Resolving this range of scales is a difficult task, but to attain a mathematically verified prediction, this breadth of scales must be captured.
Nonetheless, detailed kinetics have been used in previous studies; many of these studies have been restricted to one dimension in the inviscid limit. Sussman (Reference Sussman1995) examined pulsating hydrogen–air detonations in the inviscid limit at ambient pressure of
$0.421~\text{atm}$
. Eckett (Reference Eckett2001) considered a similar mixture of hydrogen–oxygen and found in the inviscid limit that detonations propagating into atmospheric pressure in one dimension needed to have a minimum of 150 points in the induction zone. Yungster & Radhakrishan (Reference Yungster and Radhakrishan2004) reported that the necessary resolution in the inviscid limit for hydrogen–air detonations in one dimension was near
$10^{-4}~\text{cm}$
. Daimon & Matsuo (Reference Daimon and Matsuo2007) studied various one-dimensional, inviscid, overdriven hydrogen–air mixture detonations and predicted that as the overdrive is lowered, the long time behaviour of the detonation becomes more complex. These results were to similar to those reported for other hydrogen-based fuel mixtures predicted by Sussman (Reference Sussman1995), Eckett (Reference Eckett2001) and Yungster & Radhakrishan (Reference Yungster and Radhakrishan2004).
Viscous effects have been included in several previous studies; the earliest studies by Hirschfelder & Curtiss (Reference Hirschfelder and Curtiss1958) and Wood (Reference Wood1963) are summarized by Fickett & Davis (Reference Fickett and Davis1979). These works focused on finding the viscous steady-state solution using a single step model in one dimension. Using the one-step kinetics model, Clarke, Kassoy & Riley (Reference Clarke, Kassoy and Riley1986), Clarke et al. (Reference Clarke, Kassoy, Meharzi, Riley and Vasantha1990), examined the effects of diffusion in one dimension on the development of a detonation from a small energy release. Gasser & Szmolyan (Reference Gasser and Szmolyan1993) were able to show the existence of steady diffusive strong detonations in the limit of weak diffusion using the one-step model. Singh et al. (Reference Singh, Rastigejev, Paolucci and Powers2001) calculated the development of an argon-diluted hydrogen–oxygen viscous detonation in one dimension after a reflected shock passes through the mixture. Lyng & Zumbrun (Reference Lyng and Zumbrun2004) studied the stability of a one-step detonation in the weak diffusion limit. Using a simplified kinetics model, Kivotides (Reference Kivotides2007) demonstrated that in a micro-channel the transverse wave behaviour of a detonation could be altered. Texier & Zumbrun (Reference Texier and Zumbrun2011), expanding on the earlier work of Lyng & Zumbrun (Reference Lyng and Zumbrun2004), demonstrated that in the weak diffusion limit a one-step detonation will undergo a transition through a Hopf bifurcation as the overdrive is lowered. While studying flame acceleration to a detonation in a channel, it has been shown that there is qualitative agreement between experiments and computations for a hydrogen–oxygen mixture (Liberman et al. Reference Liberman, Ivanov, Kiverin, Kuznetsov, Chukalovsky and Rakhimova2010; Ivanov, Kiverin & Liberman Reference Ivanov, Kiverin and Liberman2011; Liberman, Kiverin & Ivanov Reference Liberman, Kiverin and Ivanov2012; Ivanov et al. Reference Ivanov, Kiverin, Yakovenko and Liberman2013). Ziegler et al. (Reference Ziegler, Deiterding, Shepherd and Pullin2011) reported that for a viscous double Mach reflection detonation of hydrogen–air a near-micrometre grid-size was not quite sufficient to resolve all of the strong shocks present. Chinnayya, Hadjadj & Ngomo (Reference Chinnayya, Hadjadj and Ngomo2013) examined the viscous behaviour of a one-step detonation in a narrow channel at low pressure after formation and showed in narrower channels the addition of viscosity can completely damp the transverse wave behaviour. Recently, Lv & Ihme (Reference Lv and Ihme2014) studied several test problems for a discontinuous Galerkin method including several detailed kinetics detonations with the inclusion of diffusion.
The goal of this paper is to examine how the long time behaviour in both the time and frequency domains evolves as the supporting piston velocity is varied; in addition, it will address how the addition of physical mass, momentum and energy diffusion affects the long time behaviour in one-dimensional detonations of mixtures modelled by detailed kinetics and multicomponent transport. The plan of the paper is as follows. In § 2, the mathematical model is presented followed by a description of the computational method and a brief description of the physical problem. This is followed by a brief validation of the model and verification of the computational method. Next, the model is used to predict the dynamics of several piston-driven hydrogen–air flows. It will be demonstrated that as the supporting piston velocity is lowered, the long time behaviour becomes more complex in a similar manner to previous inviscid studies of hydrogen-based overdriven detonations (Sussman Reference Sussman1995; Eckett Reference Eckett2001; Yungster & Radhakrishan Reference Yungster and Radhakrishan2004; Daimon & Matsuo Reference Daimon and Matsuo2007). Moreover, the long time behaviour is examined using harmonic analysis and with the use of the fine resolutions in this study, gives further insight into how the behaviour changes as the supporting piston velocity is varied. In addition, this analysis is used to find similarities and differences between the inviscid and viscous analogues.
2. Formulation
2.1. Mathematical model
The governing equations are the one-dimensional, reactive, compressible Navier–Stokes equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}}{\partial t}+\frac{\partial }{\partial x}\left({\it\rho}u\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left({\it\rho}u\right)+\frac{\partial }{\partial x}\left({\it\rho}u^{2}+p-{\it\tau}\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left({\it\rho}\left(e+\frac{u^{2}}{2}\right)\right)+\frac{\partial }{\partial x}\left({\it\rho}u\left(e+\frac{u^{2}}{2}\right)+\left(p-{\it\tau}\right)u+q\right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial t}\left({\it\rho}Y_{i}\right)+\frac{\partial }{\partial x}\left({\it\rho}uY_{i}+j_{i}\right)=\overline{M}_{i}\dot{{\it\omega}}_{i},\quad i=1,\dots ,N-1, & \displaystyle\end{eqnarray}$$
$t$
, time, and
$x$
, the spatial coordinate. The dependent variables are the mixture mass density,
${\it\rho}$
, the mixture particle velocity,
$u$
, the mixture pressure,
$p$
, the viscous stress,
${\it\tau}$
, the specific internal energy of the mixture,
$e$
, the total heat flux,
$q$
, the mass fraction,
$Y_{i}$
, the diffusive mass flux,
$j_{i}$
, the molecular mass,
$\overline{M}_{i}$
and the molar production rate per unit volume,
$\dot{{\it\omega}}_{i}$
, for the
$i$
th specie.
To close the system, constitutive relations must be specified. The following constitutive relations have been chosen for an ideal mixture of calorically imperfect gases consisting of
$N$
species composed of
$L$
elements interacting in
$J$
reactions,
$$\begin{eqnarray}\displaystyle & \displaystyle j_{i}=-\frac{D_{i}^{T}}{T}\frac{\partial T}{\partial x}+{\it\rho}\mathop{\sum }_{\begin{array}{@{}c@{}}k=1\\ k\neq i\end{array}}^{N}\frac{\overline{M}_{i}D_{ik}Y_{k}}{\overline{M}}\left(\frac{1}{y_{k}}\frac{\partial y_{k}}{\partial x}+\left(1-\frac{\overline{M}_{k}}{\overline{M}}\right)\frac{1}{p}\frac{\partial p}{\partial x}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}=\frac{4}{3}{\it\mu}\frac{\partial u}{\partial x}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle q=-\mathsf{k}\frac{\partial T}{\partial x}+\mathop{\sum }_{i=1}^{N}j_{i}h_{i}-\mathscr{R}T\mathop{\sum }_{i=1}^{N}\frac{D_{i}^{T}}{\overline{M}_{i}}\left(\frac{1}{y_{i}}\frac{\partial y_{i}}{\partial x}+\left(1-\frac{\overline{M}_{i}}{\overline{M}}\right)\frac{1}{p}\frac{\partial p}{\partial x}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle y_{i}=\frac{\overline{M}}{\overline{M}_{i}}Y_{i}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \overline{M}=\left(\mathop{\sum }_{i=1}^{N}\frac{Y_{i}}{\overline{M}_{i}}\right)^{-1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle p={\it\rho}\frac{\mathscr{R}}{\overline{M}}T, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h_{i}=h_{i}^{ref}+\int _{T_{ref}}^{T}c_{P,i}(\hat{T})\text{d}\hat{T}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle e=\mathop{\sum }_{i=1}^{N}Y_{i}\left(h_{i}-\frac{\mathscr{R}T}{\overline{M}_{i}}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \dot{{\it\omega}}_{i}=\mathop{\sum }_{j=1}^{J}{\it\nu}_{ij}r_{j}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle r_{j}=k_{j}\left(\mathop{\prod }_{i=1}^{N}\left({\it\rho}\frac{Y_{i}}{\overline{M_{i}}}\right)^{{\it\nu}_{ij}^{\prime }}-\frac{1}{K_{j}^{c}}\mathop{\prod }_{i=1}^{N}\left({\it\rho}\frac{Y_{i}}{\overline{M_{i}}}\right)^{{\it\nu}_{ij}^{\prime \prime }}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle k_{j}=a_{j}T^{{\it\beta}_{j}}\exp \!\left(\frac{-\overline{E}_{j}}{\mathscr{R}T}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle K_{j}^{c}=\left(\frac{p^{ref}}{\mathscr{R}T}\right)^{\mathop{\sum }_{i=1}^{N}{\it\nu}_{ij}}\exp \!\left(\frac{-\displaystyle \mathop{\sum }_{i=1}^{N}\overline{g}_{i}^{o}{\it\nu}_{ij}}{\mathscr{R}T}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{i=1}^{N}Y_{i}=1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{i=1}^{N}j_{i}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{i=1}^{N}{\it\phi}_{li}{\it\nu}_{ij}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\nu}_{ij}={\it\nu}_{ij}^{\prime \prime }-{\it\nu}_{ij}^{\prime }, & \displaystyle\end{eqnarray}$$
$\overline{M}$
is the mixture molecular mass,
$D_{ik}$
the multi-component diffusion coefficient between the
$i$
th and
$k$
th species,
$y_{i}$
the mole fraction of the
$i$
th species,
$D_{i}^{T}$
the thermal diffusion coefficient of the
$i$
th species,
$T$
the temperature,
${\it\mu}$
the dynamic viscosity of the mixture,
$\mathsf{k}$
the thermal conductivity of the mixture,
$c_{P,i}$
the specific heat at constant pressure,
$h_{i}$
the specific enthalpy, and
$h_{i}^{ref}$
the specific enthalpy evaluated at the reference pressure of the
$i$
th species,
$\mathscr{R}$
the universal gas constant which is
$8.314\times 10^{7}~\text{erg}~\text{mol}^{-1}~\text{K}^{-1}$
,
$p^{ref}$
the reference pressure which is
$1.01325\times 10^{6}~\text{dyn}~\text{cm}^{-2}$
,
$T^{ref}$
the reference temperature which is 298.15 K,
${\it\nu}_{ij}$
the net stoichiometric coefficient,
${\it\nu}_{ij}^{\prime }$
the forward stoichiometric coefficient and
${\it\nu}_{ij}^{\prime \prime }$
the reverse stoichiometric coefficient of the
$i$
th species in the
$j$
th reaction,
$r_{j}$
the reaction rate,
$k_{j}$
the Arrhenius model for the reactions’ temperature sensitivity,
$a_{j}$
the collision frequency factor,
${\it\beta}_{j}$
the temperature exponent,
$\overline{E}_{j}$
the activation energy and
$K_{j}^{c}$
the equilibrium constant for the
$j$
th reaction,
$\overline{g}_{i}^{o}$
the
$i$
th species’ chemical potential at the reference pressure and
${\it\phi}_{li}$
the number of atoms of element
$l$
in the
$i$
th species. Equations (2.5)–(2.7) give the constitutive relations for mass, momentum and energy diffusion, which are an extended Fick’s law, a Newtonian stress–strain rate relation under Stokes’ first assumption and an extended Fourier’s law, respectively. The extended Fick’s law consists of three driving influences, the material diffusion, the pressure diffusion and the temperature diffusion, which is commonly referred to as the Soret effect. In addition, the extended Fourier’s law also consists of three driving influences, the heat flux due to heat diffusion, mass diffusion and the DuFour effect due to the pressure and material gradients. The forms of both Fourier’s and Fick’s laws are appropriate for a mixture of ideal gases, as detailed in a derivation by Merk (Reference Merk1959) and summarized by Kee et al. (Reference Kee, Dixon-Lewis, Warnatz, Coltrin and Miller1991). The reaction and mixture properties are evaluated using the CHEMKIN (Kee, Rupley & Miller Reference Kee, Rupley and Miller1992) and TRANSPORT (Kee et al.
Reference Kee, Dixon-Lewis, Warnatz, Coltrin and Miller1991) packages, respectively.
To initiate a detonation in an initially quiescent fluid, an accelerating piston is positioned on the left side of the domain. For computational purposes it is easier to use a domain attached to the face of the accelerating piston. Assuming the piston is initially located at
$x=0$
and the velocity of the piston is a known function of time,
$u_{p}(t)$
, the accelerating frame can be related to the laboratory frame as
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{x}=x-\int _{0}^{t}u_{p}(\hat{t})\text{d}\hat{t}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \tilde{t}=t, & \displaystyle\end{eqnarray}$$
$\tilde{x}$
is the location in the accelerating frame,
$\tilde{t}$
the time in the accelerating frame of reference and
$\hat{t}$
a dummy variable. Thus, the velocity in the accelerating frame,
$\tilde{u}$
, can be related to the laboratory frame particle velocity as 
$$\begin{eqnarray}\tilde{u} =u-u_{p}(t).\end{eqnarray}$$
Applying this transformation, the governing equations (2.1)–(2.4) in the accelerating frame of reference become
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\rho}}{\partial \tilde{t}}+\frac{\partial }{\partial \tilde{x}}\left({\it\rho}\tilde{u} \right)=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial \tilde{t}}\left({\it\rho}\tilde{u} \right)+\frac{\partial }{\partial \tilde{x}}\left({\it\rho}\tilde{u} ^{2}+p-{\it\tau}\right)=-{\it\rho}\frac{\text{d}u_{p}}{\text{d}\tilde{t}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial \tilde{t}}\left({\it\rho}\left(e+\frac{\tilde{u} ^{2}}{2}\right)\right)+\frac{\partial }{\partial \tilde{x}}\left({\it\rho}\tilde{u} \left(e+\frac{\tilde{u} ^{2}}{2}\right)+\left(p-{\it\tau}\right)\tilde{u} +q\right)=-{\it\rho}\tilde{u} \frac{\text{d}u_{p}}{\text{d}\tilde{t}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial \tilde{t}}\left({\it\rho}Y_{i}\right)+\frac{\partial }{\partial \tilde{x}}\left({\it\rho}\tilde{u} Y_{i}+j_{i}\right)=\overline{M}_{i}\dot{{\it\omega}}_{i},\quad i=1,\dots ,N-1. & \displaystyle\end{eqnarray}$$
$t$
will be used in place of
$\tilde{t}$
.
2.2. Computational method
The viscous calculations are performed using the wavelet adaptive multiresolution representation (WAMR) method, first developed by Vasilyev & Paolucci (Reference Vasilyev and Paolucci1996, Reference Vasilyev and Paolucci1997). This adaptive mesh refinement technique is enabled by wavelet functions. These functions have compact support in scale and space, allowing for a large compression of data. Therefore, to accurately represent a flow field, many fewer points are needed compared with a wide variety of other approaches. The WAMR method is a method of lines approach at collocation points and utilizes central finite difference schemes to calculate derivatives with special one-sided differences near boundaries. In addition, it utilizes a user-specified control threshold parameter that correlates to the error tolerated in the solution, allowing unnecessary points to be discarded. See Paolucci, Zikoski & Wirasaet (Reference Paolucci, Zikoski and Wirasaet2014b ) for a more detailed description of the method in its current form. It has been applied successfully to a number of fluids problems, see Paolucci, Zikoski & Grenga (Reference Paolucci, Zikoski and Grenga2014a ) and the references therein. The temporal integration is accomplished using an error-controlled nominally fifth-order Runge–Kutta scheme (Press et al. Reference Press, Teukolsky, Vetterling and Flannery1996).
In addition, inviscid calculations are compared directly with the viscous calculations. For these inviscid calculations, all diffusion coefficients, viscosity and thermal conductivity are taken to be zero. A uniform finite difference grid is used for these calculations and utilizes a combination of a nominally second-order mid–mod and Lax–Friedrichs scheme to calculate derivatives in a similar manner to that Xu, Aslam & Stewart (Reference Xu, Aslam and Stewart1997) implemented in their WENO and Lax–Friedrichs scheme. Temporal integration for the inviscid calculations is accomplished using a third-order Runge–Kutta method.
2.3. Physical parameters
In this study, a series of one-dimensional, piston-driven flows of an initially stoichiometric mixture of hydrogen–air (
$2\text{H}_{2}+\text{O}_{2}+3.76\text{N}_{2}$
) at ambient conditions of 293.15 K and 1 atm are considered. The detailed kinetics mechanism employed is drawn from Miller et al. (Reference Miller, Mitchell, Smooke and Kee1982) and is used by Powers & Paolucci (Reference Powers and Paolucci2005). It contains 9 species, 3 elements and 19 reversible reactions where nitrogen is treated as a non-reacting species and is shown in table 1.
The supporting piston velocity versus time curve for
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
Hydrogen–air reaction mechanism.
The flow is accelerated using a piston with a velocity,
$u_{p}(t)$
, taking the form
$$\begin{eqnarray}u_{p}(t)={\textstyle \frac{1}{2}}\left(u_{p_{o}}\left(1+\tanh \left(a\left(t-t_{a}\right)\right)\right)-\left(u_{p_{o}}-\overline{u_{p}}\right)\left(1+\tanh \left(b\left(t-t_{b}\right)\right)\right)\right),\end{eqnarray}$$
where
$u_{p_{o}}$
is the initial plateau in piston velocity,
$\overline{u_{p}}$
the final piston velocity,
$a$
the rate of acceleration to the initial plateau,
$t_{a}$
the time at which the acceleration is centred,
$b$
the rate of acceleration to the final supporting piston velocity and
$t_{b}$
the time at which the deceleration is centred. This form has two plateaus and is chosen to allow for a more rapid formation of the initial detonation. Figure 1 shows the chosen form of the piston velocity for case A in table 2 for
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
The detonation pressure versus time curve
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
. (Case A, solid light grey line; Case B, dashed black line; and Case C, dashed dark grey line.)
Initialization parameters.
The form of the piston acceleration in the initialization of the detonation plays a significant role in the early time behaviour, also known as the deflagration-to-detonation transition problem. The compression wave pushed into the fluid by the piston causes a shock to form which then proceeds to propagate away from the piston face. The strength of this shock wave is dependent on the supporting piston velocity and gives rise to a reaction wave due to preheating. This reaction wave eventually gives rise to localized explosion, which eventually develops into an overdriven detonation before relaxing at long times. This process is similar to that described by a localized thermal power deposition used by Kassoy et al. (Reference Kassoy, Kuehn, Nabity and Clarke2008), Kassoy (Reference Kassoy2010) and Regele, Kassoy & Vasilyev (Reference Regele, Kassoy and Vasilyev2012). The weaker the inertial confinement of the initial reaction wave, the longer the initial detonation takes to form.
Because this work focuses on the long time behaviour of the detonation, the long time behaviour of three piston acceleration profiles is examined with the same final supporting piston velocity of
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
to examine whether this early time effect continues to later times. The parameters for the three different initialization profiles are listed in table 2. As shown in figure 2, the inertial confinement of the initial reaction wave plays a dramatic role in the time to detonation. In addition, the over-pressure is much greater for the weaker driven shock, and it also takes longer to relax to a steady state. However, the long time behaviour relaxes to a stable detonation in all three cases, with the local maximum pressure at the front, which is the detonation pressure, of
$p_{A}=36.75~\text{atm}$
,
$p_{B}=36.76~\text{atm}$
and
$p_{C}=36.73~\text{atm}$
, meaning all three cases are within 0.1 % of each other. For the remainder of the paper, case A is used allowing only
$\overline{u_{p}}$
to vary. Final piston velocities ranging from
$1.200\times 10^{5}$
to
$1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
are examined.
For the viscous calculations the threshold used for the WAMR method is
${\it\epsilon}=10^{-3}$
unless otherwise specified. This selection leads to a spatial resolution of
$O(3\times 10^{-6}~\text{cm})$
to be utilized which results in a time step of
$O(10^{-12}~\text{s})$
. After the formation of the detonation, a typical simulation time for a viscous calculation of
${\sim}1~{\rm\mu}\text{s}$
took on
$O(300~\text{CPUh})$
on 32 cores. For the inviscid calculations a spatial discretization of
$2.5\times 10^{-5}~\text{cm}$
is used which results in a time-step of
$O(10^{-10}~\text{s})$
.
3. Validation and verification
As the model has been restricted to one dimension, there are limited means of validation. However, in experiments of shock-induced combustion flow around spherical projectiles in a hydrogen–air mixture at 0.421 atm and 293.15 K, Lehr (Reference Lehr1972) observed longitudinal oscillations. For an inflow condition corresponding to an overdrive
$f\approx 1.10$
, Lehr observed a frequency of
${\it\nu}_{o}=1.04~\text{MHz}$
, where the overdrive
$f=(D/D_{CJ})^{2}$
and
$D$
and
$D_{CJ}$
are the detonation wave speed and the wave speed at which the detonation terminates at the sonic point, respectively. Starting the viscous calculation with the inviscid steady-state profile with a superimposed smooth transition from the shocked state to the ambient condition over
$5\times 10^{-4}~\text{cm}$
, a frequency of
${\it\nu}_{o}=0.97~\text{MHz}$
is predicted. Thus, it seems that the instability observed by Lehr in multiple dimensions is captured well by a one-dimensional model. This is similar to results reported by Yungster & Radhakrishan (Reference Yungster and Radhakrishan2004) of
${\it\nu}_{o}=1.06~\text{MHz}$
for an overdrive of
$f=1.09$
with an ambient temperature of 298 K and Daimon & Matsuo (Reference Daimon and Matsuo2007) who do not report frequency explicitly, but visual inspection suggests that the frequencies are in the 1 MHz range. The predicted frequency here is only 6.7 % different from that measured by Lehr; the discrepancy is likely due the one-dimensional assumption and uncertainty in chemical kinetic parameters.
The WAMR procedure is a self-converging method, which means that as the error-threshold parameter,
${\it\epsilon}$
, is reduced, the overall error is reduced. To verify that in fact the procedure is convergent regime, several values of
${\it\epsilon}$
are examined for
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
. The long time behaviour at this supporting piston velocity is a stable, steadily travelling detonation as shown in figure 2. Figure 3(a) shows the long time detonation pressure versus
${\it\epsilon}$
. The detonation pressure is converging to 36.68 atm. It should be noted that the standard deviation on the detonation pressure are indicated by the vertical lines. As the error-threshold parameter is reduced, the standard deviation is reduced around the detonation pressure point indicated by the dots. In fact, at the two most accurate solutions, the standard deviation in the detonation pressure is difficult to identify. In addition, the difference in the long time detonation pressure is calculated from the most accurate solution; this is shown in figure 3(b). As the error-threshold parameter is reduced the difference decreases near
$O({\it\epsilon}^{0.9})$
, as indicated on the log–log plot. Furthermore, the largest percentage difference is 0.2 % giving a good indication that the WAMR method is in the convergent regime.
(a) Detonation pressure versus
${\it\epsilon}$
and (b) difference in detonation pressure between
${\it\epsilon}=1\times 10^{-6}$
and
${\it\epsilon}$
for
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
4. Results
In this section, a study of the long time behaviour of the propagating detonation is performed as the final supporting piston velocity,
$\overline{u_{p}}$
, is varied. This is done first in the time domain, and then harmonic analysis is used to examine the active frequencies of the pulsating detonations. Several comparisons between the viscous and inviscid calculations are performed in both the time and frequency domains.
(a) Several time shots of the spatial pressure profile (solid black line,
$10\times 10^{-6}~\text{s}$
; solid light grey line,
$35\times 10^{-6}~\text{s}$
; and dashed grey line,
$60\times 10^{-6}~\text{s}$
) and (b) typical spatial profile of mass fractions at a
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
4.1. Stable detonations
For sufficiently high
$\overline{u_{p}}$
, a steadily travelling detonation arises and persists at long times. The detonation pressure versus time curve for a stable detonation, at
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
, is shown in figure 1. For case A, by
$10\times 10^{-6}~\text{s}$
the detonation relaxes to a steadily travelling piston-supported detonation travelling to the right at
$2.244\times 10^{5}~\text{cm}~\text{s}^{-1}$
. Spatial pressure profiles after the detonation relaxes to the stable detonation are shown in figure 4(a) for
$t=10\times 10^{-6}~\text{s}$
,
$t=35\times 10^{-6}~\text{s}$
and
$t=60\times 10^{-6}~\text{s}$
. The later time profiles have been shifted in space using the steady wave speed. There are only minuscule differences between the front locations; these differences are more clearly shown in the inset. However, the largest difference between the front locations is still only
$2.5\times 10^{-4}~\text{cm}$
. Figure 4(b) shows the spatial mass fraction at
$t=50\times 10^{-6}~\text{s}$
which is representative of the steadily travelling detonation front. As a particle passes through the detonation, it first encounters a thin viscous shock accompanied by rapid pressure and temperature rise. Then, its pressure and temperature remain relatively constant as it traverses a short induction zone. In this zone radicals are generated. When a sufficient number of radicals are present, the fluid particle enters a thin zone in which vigorous reaction commences. Here pressure and temperature vary rapidly. Finally, it passes into a thick relaxation zone, where all state variables equilibrate.
4.2. High-frequency pulsating detonations
After
$\overline{u_{p}}$
is lowered below a critical value, the long time behaviour of the propagating detonation undergoes a transition from a steadily travelling wave to a pulsating detonation. This transition occurs between
$\overline{u_{p}}=1.420\times 10^{5}~\text{cm}~\text{s}^{-1}$
and
$\overline{u_{p}}=1.410\times 10^{5}~\text{cm}~\text{s}^{-1}$
. Figure 5(a) shows the detonation pressure versus time curves for a supporting piston velocity just above and just below the transition point. These pulsations are caused by the slight detachment between the pressure wave and the reaction wave, which in turn elongates the induction zone. The phase space plot for both the stable and unstable case for both the stable and unstable case is shown in figure 5(b). For the stable case, the phase space plot is a black dot located at
$p=34.95~\text{atm}$
, the dot is enlarged by 10 times for ease of viewing. At
$\overline{u_{p}}=1.410\times 10^{5}~\text{cm}~\text{s}^{-1}$
, it becomes clear that detonation is pulsating; however, it is difficult to extract whether it is near cyclic from the phase space plot. This case will be examined further in § 4.5 to extract more information about the pulsations.
(a) Detonation pressure versus time and (b) phase space plot for both
$\overline{u_{p}}=1.420\times 10^{5}~\text{cm}~\text{s}^{-1}$
(dashed black line) and
$\overline{u_{p}}=1.410\times 10^{5}~\text{cm}~\text{s}^{-1}$
(solid grey line). Note the phase space plot for
$\overline{u_{p}}=1.420\times 10^{5}~\text{cm}~\text{s}^{-1}$
has been enlarged by 10 times.
In contrast to the clear periodic limit cycles predicted by Henrick, Aslam & Powers (Reference Henrick, Aslam and Powers2006) for the simple one-step model in the CJ limit, the pulsating detonations here do not produce nearly as smooth limit cycles. This is likely influenced by several factors. First, the piston-driven flows in this study are overdriven in nature; as such, the positively moving characteristic waves travel through different decaying
$N$
-waves in the negatively moving characteristic field emanating from the detonation front. The likelihood of these positively moving characteristic waves and decaying negatively moving characteristic
$N$
-waves being synchronized is extremely low, and thus precludes precisely periodic cycles. These positively moving characteristics in the overdriven case clearly reach the detonation shock front. In the CJ case, there is a sonic locus that remains a finite distance behind the front. As demonstrated by Kasimov & Stewart (Reference Kasimov and Stewart2004) for the one-step model, this sonic locus acts as an information barrier. It only allows characteristics in front of it to propagate towards the front. In addition, the one-step model has only a single length scale of reaction, whereas the detailed hydrogen–air mechanism has reaction length scales that span several orders of magnitude. Furthermore, the one-step model is irreversible, while the detailed kinetics has reversible reactions.
The induction zone length changes in the pulsating detonations due to slight separation between the pressure wave and the reaction front. For
$\overline{u_{p}}=1.410\times 10^{5}~\text{cm}~\text{s}^{-1}$
, the minimum induction zone length is
$8.06\times 10^{-3}~\text{cm}$
during the cycle which corresponds to the point of peak detonation pressure. The maximum temperature gradient behind the front is used as the indicator for the end of the induction zone length. The maximum induction zone length occurs at the minimum detonation pressure in the cycle and is
$8.73\times 10^{-3}~\text{cm}$
. As figure 6 shows, the stable detonation just above the bifurcation point, at
$\overline{u_{p}}=1.420\times 10^{5}~\text{cm}~\text{s}^{-1}$
, has a induction zone length of
$8.11\times 10^{-3}~\text{cm}$
. This induction zone length is close to the minimum of that of the pulsating detonation because the detonation pressure for the stable case is similar to the peak detonation pressure in the pulsating case. At the minimum detonation pressure in the cycle, it is clear that the gradient in temperature is delayed and decreased in magnitude; likewise, at the peak detonation pressure the peak is slightly greater than that in the stable case.
The temperature gradient through the induction zone for a stable detonation at
$\overline{u_{p}}=1.420\times 10^{5}~\text{cm}~\text{s}^{-1}$
(black solid line), at the peak detonation pressure (dashed light grey line) and at the minimum detonation pressure (dark grey line) for an unstable detonation at
$\overline{u_{p}}=1.410\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
As shown in figure 7(a), when
$\overline{u_{p}}$
is lowered further below the bifurcation point, the oscillations grow in amplitude. In addition, the frequency shifts towards lower frequencies. As the pulsations become larger in amplitude, it becomes clearer that they are nearly periodic as demonstrated by successive pulsations nearly coinciding in the phase plot of figure 7(b).
(a) Detonation pressure versus time and (b) phase space plot for both
$\overline{u_{p}}=1.400\times 10^{5}~\text{cm}~\text{s}^{-1}$
(dashed black line) and
$\overline{u_{p}}=1.370\times 10^{5}~\text{cm}~\text{s}^{-1}$
(solid grey line).
4.3. Multiple mode pulsating detonations
The behaviour becomes even more complex at lower supporting piston velocities with a dual mode behaviour arising below a second bifurcation point. An example of this type of propagating detonation is shown in figure 8(a) for
$\overline{u_{p}}=1.310\times 10^{5}~\text{cm}~\text{s}^{-1}$
. It is apparent that the dual mode behaviour persists at long times. Although these dual mode detonations do not repeat in a clean limit cycle, it is still obvious that it is stably bounded at long times which is demonstrated by the phase plane plot shown in figure 8(b).
(a) Detonation pressure versus time and (b) phase space plot for a
$\overline{u_{p}}=1.310\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
4.4. Low-frequency-dominated pulsating detonations and chaos
At yet even lower supporting piston velocities, this dual mode behaviour relaxes into a mode that is dominated by a single low-frequency pulsating flow. Figure 9(a) shows this relaxation to a nearly periodic limit cycle at long times. However, the phase space plot, shown in figure 9(b), indicates that even at long times there is still some variation in the cycle. Once this low-frequency mode becomes the dominant mode, a behaviour similar to period-doubling is predicted. As shown in figure 9(c), a nearly period-two detonation is predicted at
$\overline{u_{p}}=1.230\times 10^{5}~\text{cm}~\text{s}^{-1}$
. At this supporting piston velocity, the relative maxima can be grouped into two distinct groups; the first at
$p_{1}=47.56\pm 0.68~\text{atm}$
and the second being
$p_{2}=50.9\pm 0.84~\text{atm}$
. Figure 9(d) clearly exhibits the distinct two-lobe phase space for a period-two detonation. This phenomenon is exhibited even more prominently at
$\overline{u_{p}}=1.220\times 10^{5}~\text{cm}~\text{s}^{-1}$
, as shown in figure 9(e). However, the higher relative maxima is more erratic as indicated by the wider spread in right-most lobe shown in the phase space plot of figure 9(f). This period-doubling behaviour is more clearly seen in the frequency domain and will be discussed further in § 4.5. After this period-doubling regime, the detonation pressure versus time curve exhibits many more relative maxima, which is shown in figure 9(g) for
$\overline{u_{p}}=1.200\times 10^{5}~\text{cm}~\text{s}^{-1}$
. This is further elucidated by examining the phase space plot, shown in figure 9(h), where no consistent cycle is visible. The system likely underwent a transition to chaos. However, to definitively categorize the system as chaotic further analysis would be needed.
Detonation pressure versus time and phase plots for (a,b)
$\overline{u_{p}}=1.250\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (c,d)
$\overline{u_{p}}=1.230\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (e, f)
$\overline{u_{p}}=1.220\times 10^{5}~\text{cm}~\text{s}^{-1}$
and (g,h)
$\overline{u_{p}}=1.200\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
PSD viscous spectra at (a)
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (b)
$\overline{u_{p}}=1.410\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (c)
$\overline{u_{p}}=1.340\times 10^{5}~\text{cm}~\text{s}^{-1}$
and (d)
$\overline{u_{p}}=1.310\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
4.5. Harmonic analysis
Next, the detonation pressure versus time behaviour is examined using harmonic analysis. This is a tool that can be used to extract important information of a signal, such as the dominant frequency and a ratio of energy carried at various frequencies. The specific tool that is used is the power spectral density (PSD). This tool is chosen because it describes how the variance (or power) is distributed in the frequency domain. It is real valued for any real signal. The PSD is defined as the Fourier transform of the auto-correlation of a signal (Hamilton Reference Hamilton1994; Billinger Reference Billinger2001) and can be written as the magnitude squared of the Fourier transform of the signal by using the Wiener–Khinchin theorem. The PSD of a signal reveals periodicities that can be hidden in a complex signal, e.g. the dual mode pulsating behaviour. Furthermore, it can be helpful to discern how the frequency of the pulsations is affected by changing the supporting piston velocity. The detonation pressure time series curves are analysed after the initialization period.
In this study, the discrete one-sided mean-squared amplitude PSD is used. The single-sided PSD is chosen so that the aliasing effect at high frequencies could be bypassed. This normalization is chosen such that, as Parseval’s theorem states (Oppenheim & Schafer Reference Oppenheim and Schafer1975), the sum of
${\it\Phi}_{d}$
to equal the mean-squared amplitude of the discrete detonation pressure signal, where
${\it\Phi}_{d}({\it\nu}_{k})$
is the discrete PSD of the detonation pressure–time signal at frequency,
${\it\nu}_{k}$
. Ng et al. (Reference Ng, Higgins, Kiyanda, Radulescu, Lee, Bates and Nikiforakis2005) applied a similar type of procedure to pulsating detonations using the simple one step kinetics model in the inviscid limit. The PSDs presented are for the detonation pressure–time signal. For signals with deviations larger than
$0.04~\text{atm}$
from the mean, the mean detonation pressure is subtracted for calculation.
PSD viscous spectra at (a)
$\overline{u_{p}}=1.290\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (b)
$\overline{u_{p}}=1.260\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (c)
$\overline{u_{p}}=1.230\times 10^{5}~\text{cm}~\text{s}^{-1}$
and (d)
$\overline{u_{p}}=1.220\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
As discussed previously in § 4.1, when the supporting piston velocity is sufficiently high, the long time behaviour is a steadily travelling detonation wave. As depicted in the frequency domain as shown in figure 10(a), the PSD spectrum demonstrates all of the energy is concentrated near the zero frequency. Lowering
$\overline{u_{p}}$
below the first bifurcation point, gives rises to a pulsation at
${\it\nu}_{o}=3.41~\text{MHz}$
at
$\overline{u_{p}}=1.410\times 10^{5}~\text{cm}~\text{s}^{-1}$
. In figure 10(b) it is clear that the majority of the pulsation energy is carried at a single frequency. However, the second harmonic frequency also carries energy; this results in slight differences in the relative maxima in detonation pressure in cycle. As the supporting piston velocity is lowered, the frequency spectrum blue shifts. At
$\overline{u_{p}}=1.340\times 10^{5}~\text{cm}~\text{s}^{-1}$
, shown in figure 10(c), the fundamental frequency is now located at
${\it\nu}_{o}=2.85~\text{MHz}$
, and the harmonics have shifted as well. In fact, the ratio of the amount of energy being concentrated at higher harmonics has increased, which is demonstrated by the appearance of the third harmonic in the plot. Examining a
$\overline{u_{p}}$
further below the neutral stability boundary, it becomes clear there is a low-frequency mode that is now playing an important role in the long time behaviour of the pulsations as shown for
$\overline{u_{p}}=1.310\times 10^{5}~\text{cm}~\text{s}^{-1}$
in figure 10(d). At this
$\overline{u_{p}}$
, the low-frequency mode occurs at
${\it\nu}_{l}=0.44~\text{MHz}$
and carries a significant amount of energy. However, the high-frequency mode, which occurs at
${\it\nu}_{h}=2.65~\text{MHz}$
, is still the dominant mode. These modes remain stationary at a smaller error-threshold parameter, once again confirming the WAMR is capturing the long time behaviour. Furthermore in this regime where there are two dominant modes; the modes interact giving rise to many more modes that carry energy.
This interaction of the two dominant modes gives rise to a modulation instability. This modulation instability phenomenon occurs in many other physical systems due to the inherent nonlinearity of the physical world (Zakharov & Ostrovsky Reference Zakharov and Ostrovsky2009). In these pulsating detonations it manifests itself as active modes, called sidebands, at multiples of the low frequency around the high frequency and its harmonics. These active modes surrounding the high-frequency mode and its harmonics form envelope waves that persist at long times.
After the appearance of the dual mode behaviour, the ratio of the energy present in the pulsation carried at the high fundamental frequency continues to decrease as
$\overline{u_{p}}$
is lowered further. Figure 11(a) shows the PSD for a
$\overline{u_{p}}=1.290\times 10^{5}~\text{cm}~\text{s}^{-1}$
and demonstrates that a much more dominant low fundamental frequency at 0.41 MHz exists compared with that of figure 10(d). However, the high-frequency mode, which is located at 2.64 MHz, still carries energy. In fact, the lower side bands around the high-frequency mode carries a similar order of magnitude of energy as the second harmonic of the low-frequency mode. However, it is also clear the side bands have also been reduced, indicating further that the high-frequency modes have weakened. Even as more energy is shifting to the lower frequency, the spectrum continues to blue shift towards lower frequencies as the supporting piston velocity is lowered, but at a slower rate than predicted in the high-frequency mode. Eventually the low frequencies become so dominant that the high-frequency mode and side bands carry less energy than the fourth harmonic of the low-frequency mode, as shown figure 11(b) by the PSD at
$\overline{u_{p}}=1.260\times 10^{5}~\text{cm}~\text{s}^{-1}$
, where
${\it\nu}_{o}=0.38~\text{MHz}$
. Nonetheless, there is a side band frequency mode that still persists at long times, but at a lower energy state. This is likely a manifestation of the multiple reaction length scales interacting with each other as well as the diffusion length scales.
As briefly mentioned in § 4.4, after the low-frequency mode has become dominant, the long time behaviour goes through a phenomenon similar to period-doubling. This is more clearly illustrated in the frequency domain shown in the two PSD spectra in figure 11(c,d) for
$\overline{u_{p}}=1.230\times 10^{5}~\text{cm}~\text{s}^{-1}$
and
$\overline{u_{p}}=1.220\times 10^{5}~\text{cm}~\text{s}^{-1}$
, respectively. This near period-doubling is illustrated by the appearance of subharmonics of the fundamental frequency. As an example, in figure 10(c) the fundamental frequency is located at 0.34 MHz, but there are peaks in the PSD spectrum at 0.17 MHz and 0.49 MHz, which
$1/2$
and
$3/2$
the fundamental frequency, respectively. These are subharmonic frequencies, which indicates that the long time behaviour of the pulsations is near a limit cycle with two distinct relative maxima in the detonation pressure time curve. Figure 11(d) shows that the first set of subharmonics have grown in amplitude indicating the strength of the period-two detonation has grown.
4.6. Bifurcation diagram
A bifurcation diagram is constructed showing the various propagation modes. It has been created with 31 supporting piston velocities spaced at
$1.0\times 10^{3}~\text{cm}~\text{s}^{-1}$
and as such is a coarse approximation of the full diagram. Figure 12(a) shows how the maximum detonation pressure evolves versus the supporting piston velocity. Note the peak detonation pressure has been scaled by the average detonation pressure. As the peak detonation pressure varies from cycle to cycle, the standard deviation of peaks in the stable, high-frequency-dominated and low-frequency-dominated modes are indicated by vertical lines. In both the stable and high-frequency-dominated modes, the standard deviations are difficult to distinguish from the peak detonation pressure. The region in which there are two active modes is indicated by the dense number of points near
$\overline{u_{p}}=1.300\times 10^{5}~\text{cm}~\text{s}^{-1}$
; likewise, the dense region near
$\overline{u_{p}}=1.200\times 10^{5}~\text{cm}~\text{s}^{-1}$
is indicative of a detonation with many active modes, which is likely chaotic. This is more clearly understood by looking at the bifurcation plot in of active frequencies, shown in figure 12(b), in which the shade of the points indicates the magnitude with the darkest being the most dominant mode and the lightest being the weakest. In the high-frequency mode, there are three active frequencies: the fundamental frequency, the second harmonic, and the third harmonic. The blue shift of the frequency spectrum is most clearly seen in the third harmonic. In the dual mode region, it is apparent that side banding occurs near the high-frequency mode and its harmonics; however, there are still just two dominant modes. The side banding continues in the low-frequency mode, but at weaker strengths than that of the dual mode. In addition, subharmonics appear at
$1/2$
and
$3/2$
at both
$\overline{u_{p}}=1.230\times 10^{5}~\text{cm}~\text{s}^{-1}$
and
$\overline{u_{p}}=1.220\times 10^{5}~\text{cm}~\text{s}^{-1}$
. At this lower supporting piston velocity, further subharmonics appear at the half intervals as well as the previously mentioned subharmonics grow in strength. In the lowest supporting piston velocities studied, many frequencies are active indicating that it is likely that the detonation is in a chaotic regime.
Bifurcation plot of (a) the maximum detonation pressure scaled by the average detonation pressure (b) active frequencies versus supporting piston velocity where the darker shade of the point indicates larger magnitude.
Detonation pressure versus time for both viscous (black lines) and inviscid (grey lines) cases at (a)
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (b)
$\overline{u_{p}}=1.430\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (c)
$\overline{u_{p}}=1.400\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (d)
$\overline{u_{p}}=1.320\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (e)
$\overline{u_{p}}=1.310\times 10^{5}~\text{cm}~\text{s}^{-1}$
and ( f)
$\overline{u_{p}}=1.250\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
4.7. Comparison with the inviscid analogue
Several supporting piston velocities are examined in the inviscid limit to elucidate the effects of physical diffusion on a one-dimensional detonation of detailed kinetics mechanism where instabilities are manifested as pulsations. As in the viscous case, when the supporting piston velocity is sufficiently high, a stable steadily travelling detonation is formed and persists at long times. Figure 13(a) shows both the viscous and inviscid detonation pressure versus time curves at
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
. The inviscid case relaxes to a detonation pressure of
$36.68~\text{atm}$
which is less than
$0.1\,\%$
different from the viscous analogue at this piston velocity. However, at
$\overline{u_{p}}=1.430\times 10^{5}~\text{cm}~\text{s}^{-1}$
the inviscid detonation begins to pulsate with an oscillation amplitude of
${\sim}1~\text{atm}$
whereas the viscous analogue remains stable as shown in figure 13(b). This pulsation amplitude is larger than that of the viscous case at
$\overline{u_{p}}=1.410\times 10^{5}~\text{cm}~\text{s}^{-1}$
. Thus, the addition of diffusion to the model has added a slightly stabilizing effect, shifting the transition to a pulsating detonation by greater than
$1.5\,\%$
, but less than 2 % with respect to the supporting piston velocity. Figure 13(c) shows the long time behaviour at
$\overline{u_{p}}=1.400\times 10^{5}~\text{cm}~\text{s}^{-1}$
for both the inviscid and viscous cases. The relative maxima in detonation pressure are
$p=35.66\pm 0.10~\text{atm}$
and
$p=36.62\pm 0.005~\text{atm}$
, for the viscous and inviscid cases, respectively. In addition to the reduction of the maximum detonation pressure, the amplitude of oscillations has also been reduced by 40 % by the addition of viscosity. However, as the pulsations become stronger, the effect of viscosity is reduced as demonstrated in figure 13(d) for
$\overline{u_{p}}=1.320\times 10^{5}~\text{cm}~\text{s}^{-1}$
. The pulsation amplitude reduction due to diffusion is weakened to less than
$0.1\,\%$
near the transition point to the dual mode behaviour. Figure 13(e) shows the detonation pressure versus time curve for both the inviscid and viscous cases at
$\overline{u_{p}}=1.310\times 10^{5}~\text{cm}~\text{s}^{-1}$
, which is in the dual mode pulsating behaviour in both cases. It is difficult to identify differences in the time domain due to the interacting modes; the frequency domain will be discussed later. The average detonation pressure and the average maximum detonation for the viscous case are
$p=34.32~\text{atm}$
and
$p=36.2\pm 0.7~\text{atm}$
. Likewise for the inviscid case, the average detonation pressure and the average maximum detonation are
$p=34.33~\text{atm}$
and
$p=36.1\pm 0.7~\text{atm}$
. In the low-frequency-dominated mode the effect of viscosity is nearly negligible, which is demonstrated in figure 13(f) for
$\overline{u_{p}}=1.250\times 10^{5}~\text{cm}~\text{s}^{-1}$
. The local maxima are
$p=46.5\pm 0.4~\text{atm}$
and
$p=46.6\pm 0.5~\text{atm}$
, for the viscous and inviscid cases; respectively. This is a relative difference of 0.2 %, and it is clear that the maxima overlap.
PSD spectra for both viscous (black lines) and inviscid (grey lines) cases at (a)
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (b)
$\overline{u_{p}}=1.430\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (c)
$\overline{u_{p}}=1.400\times 10^{5}~\text{cm}~\text{s}^{-1}$
,(d)
$\overline{u_{p}}=1.320\times 10^{5}~\text{cm}~\text{s}^{-1}$
, (e)
$\overline{u_{p}}=1.310\times 10^{5}~\text{cm}~\text{s}^{-1}$
and (f)
$\overline{u_{p}}=1.250\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
The PSD spectra are calculated and compared using the average inviscid detonation pressure to scale both the inviscid and viscous detonation pressures; furthermore, the PSD is calculated in decibels using the maximum value of either case. Supporting piston velocities ranging from
$1.250\times 10^{5}$
to
$1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
are shown in figure 14. At
$\overline{u_{p}}=1.500\times 10^{5}~\text{cm}~\text{s}^{-1}$
, both of the PSD spectra are concentrated around the zero frequency, as shown in figure 14(a), indicating that the detonation is stable at long times. Figure 14(b) shows the PSD spectra at
$\overline{u_{p}}=1.430\times 10^{5}~\text{cm}~\text{s}^{-1}$
; it is clear the inviscid case is pulsating at
${\it\nu}_{o}=3.60~\text{MHz}$
, but the PSD of the viscous calculation is still concentrated around the zero frequency indicating a stable detonation. At
$\overline{u_{p}}=1.400\times 10^{5}~\text{cm}~\text{s}^{-1}$
, the fundamental frequency in the inviscid case is
${\it\nu}_{o}=3.36~\text{MHz}$
, whereas in the viscous case it is minimally shifted to a lower frequency by 1 %. The shift is more apparent in the second harmonic, which is shifted by
$0.06~\text{MHz}$
. Figure 14(c) shows that the magnitude of fundamental frequency is larger in the inviscid case indicating the reduction in pulsation amplitude. The addition of viscosity affects more the size of the pulsation than the frequency of pulsations. Near the transition point to the dual mode behaviour, at
$\overline{u_{p}}=1.320\times 10^{5}~\text{cm}~\text{s}^{-1}$
, the frequency shift is reduced to
$0.01~\text{MHz}$
. Figure 14(d) indicates that the fundamental frequency peaks in the inviscid and viscous case are closer in magnitude than at the higher supporting piston velocities, giving another indication that the pulsation amplitude is nearly identical. In the dual mode, as shown in figure 14(e) the active modes are only barely distinguishable from each other; however, the strength of the low-frequency mode is stronger in the inviscid case, and the high-frequency mode is weaker. This indicates that, though small, the addition of physical viscosity to the model is still playing a role and slightly delays the transition to the dual mode behaviour. In addition, the high-frequency mode is shifted, but only by 0.3 %. When a low-frequency dominated mode (
$\overline{u_{p}}=1.250\times 10^{5}~\text{cm}~\text{s}^{-1}$
) is examined, it is seen that the PSD spectra, shown in figure 14(f), are nearly indistinguishable from each other. The fundamental frequencies are identical, and as the magnitude of the pulsations are the same, the magnitude of the PSD at this frequency are also identical. However, the PSD for the viscous case is missing the fifth harmonic and has minimally more energy carried at the higher frequencies.
The amplitude reduction present in the high-frequency mode is weakened as the supporting piston velocity is lowered. At lower piston velocities the intrinsic instability grows stronger and, thus, the effect of physical viscosity is weaker. The addition of physical viscosity to the model has an overall stabilizing effect, delaying the initial transition to instability and reducing the amplitude of oscillations in the pulsating mode dominated by high-frequency oscillations. This suggests that in multiple dimensions, that diffusion can play an important role in the formation and propagation in detonations in narrow channels, where the transverse waves can possibly be damped. This will examined with detailed kinetics in a future study. However, the further away from this transition, physical viscosity plays a less important role in determining the long time behaviour at least in one dimension.
5. Discussion
It is useful to consider a physics-based interpretation of these detailed unsteady detonation dynamics. Our interpretation is supported either directly by the current results or plausible hypotheses that could guide future studies. The discussion will be mainly cast in the framework of the current one-dimensional piston-supported detonations in a viscous hydrogen–air mixture; to quantitatively illustrate these important points for viscous limit cycle detonations, an appeal is made to a simpler model of one-step kinetics where the CJ limit cycle is well quantified (Romick et al. Reference Romick, Aslam and Powers2012). It is well understood that the compressible reactive Navier–Stokes equations admit steady travelling wave solutions in response to a driving piston. The steady wave is driven by a combination of mechanical energy input from the driving piston and chemical energy input from the exothermic heat release in the subsonic region following the thin lead shock. For sufficiently high piston speeds such that the kinetic energy imparted by the piston is much greater than the chemical energy, the wave behaves similarly to an inert shock wave. As the piston velocity is lowered, the chemical energy makes an ever-increasing relative contribution to driving the wave. At a critical CJ piston velocity, all of the energy to drive the wave is available from the chemical energy, and the wave becomes self-propagating.
These steady waves can respond differently to small perturbations in the various regimes of supporting piston velocity. The key question is whether such perturbations grow or decay, and if they initially grow, what physical mechanisms prevent unbounded growth. In general terms, there are two physical mechanisms present which induce dissipation of structured mechanical and chemical energy into unstructured thermal energy: diffusion and the irreversible part of chemical reaction. Simultaneously, there are physical mechanisms present which induce the growth and resonance of various oscillatory structures predicted in some cases: amplification of selected modes by exothermic reaction combined with the effects of advection and diffusion. Ultimately, nonlinearity has a role. It is often the case that modes which grow linearly away from equilibrium can move into a region where nonlinear effects become important and serve to either suppress further growth or induce some variety of catastrophic growth. The action of these various physical mechanisms is a strong function of the various length and time scales in play as the driving piston velocity is varied.
For sufficiently high piston velocity, the effect of exothermic chemistry is minimal. One might imagine that a small sinusoidal disturbance near the shock front would segregate into one entropic and two acoustic modes, travelling near the local particle and acoustic speeds, respectively. Diffusion would act to reduce the amplitude of the disturbance. High-frequency modes would dissipate more rapidly than low-frequency modes, but ultimately all would stabilize, and the system would relax to a steady propagating wave. Such is what is predicted for
$\overline{u}_{p}>1.420\times 10^{5}~\text{cm}~\text{s}^{-1}$
.
For lower piston velocities, e.g.
$1.400\times 10^{5}~\text{cm}~\text{s}^{-1}$
, it is obvious that limit-cycle-like behaviour is predicted. Thus, at this piston velocity, nature favors a partition of the chemical and kinetic energy of the fluid into a pattern in which some of the energy resides in the two modes displayed in figure 14(c). Under these conditions, there is little difference between the viscous and inviscid predictions, so it is inferred that the physics are best understood in the context of a reactive Euler model. This piston velocity is likely favourable for the establishment of an organ-pipe-type resonance influenced by a balance between reaction and advection. The relevant length is the induction zone, that is, the region between the lead shock and the point where significant chemical reaction commences. Good estimates for a similar mixture are given by Powers & Paolucci (Reference Powers and Paolucci2005). The reaction kinetics are such that the induction zone length
$\ell _{ind}\simeq 10^{-2}~\text{cm}$
. And the material properties are such that the post-shock acoustic speed
$c\simeq 10^{5}~\text{cm}~\text{s}^{-1}$
. A rough estimate of the fundamental resonant frequency is thus
${\it\nu}\simeq c/\ell _{ind}=10~\text{MHz}$
. This is of the same order of magnitude as that predicted. This scaling argument is consistent with the results of Short (Reference Short1996), who showed that perturbations within the induction zone were linearly unstable while examining one-step square wave detonations at high activation energy. Moreover, for this case, figure 14(c) reveals that diffusion induces a small-amplitude reduction in the resonant modes, as well as a small shift in the resonant frequencies. This is consistent with what is found in ordinary nonlinear mass–spring–damper systems (Strogatz Reference Strogatz2014). Again, similar to what is found in a nonlinear mass–spring–damper system, it is most likely that nonlinear effects serve a much stronger role in suppressing the growth of the resonant modes. For higher-frequency modes, it is likely that diffusion plays the dominant role in amplitude suppression. With the finest length scale of reaction given by Powers & Paolucci (Reference Powers and Paolucci2005) for this mixture
$\ell _{finest}\simeq 10^{-4}~\text{cm}$
, and the diffusivity of the mixture
$D_{mix}\simeq 10~\text{cm}^{2}~\text{s}^{-1}$
, one can estimate the frequency of the disturbance for which diffusion clearly dominates as
${\it\nu}\simeq D_{mix}/\ell _{finest}^{2}=1~\text{GHz}$
. Useful insights on the relative importance of advection, reaction and diffusion in hydrogen–air chemistry is given by Al-Khateeb et al. (Reference Al-Khateeb, Powers and Paolucci2013) in the context of a laminar flame. There, it is shown that diffusion clearly influences the various reaction time scales on length scales given by a classical Maxwellian model,
$\ell _{i}\simeq \sqrt{D_{mix}{\it\tau}_{i}}$
, where
$\ell _{i}$
is the reaction length scale associated with the chemical time scale
${\it\tau}_{i}$
.
As the piston velocity is lowered further, nonlinearity plays a more prominent role, especially as seen in figures 10(d) and 11(a–d). As the oscillatory modes are dominated here by ever-lower frequencies as the piston velocity is lowered, it is likely that diffusion is playing even less of a role in the dominant low-frequency dynamics, with its main effect being confined to much higher-frequency modes. Even then, the presence of diffusion is important in providing a physically based cutoff mechanism for high-frequency modes. As documented by Powers (Reference Powers2006), lack of such a cutoff mechanism can then admit approximations which do not converge as the grid discretization scale is reduced, thus rendering the results to be potentially strongly influenced by the size of the discretization and the selected numerical method. Moreover, Mazaheri et al. (Reference Mazaheri, Mahmoudi and Radulescu2012) demonstrated that, in regions that have large gradients in the flow, diffusion plays an influential role.
As it is prohibitively expensive to fully relax detonations with detailed kinetics to limit cycles, these notions of piston-driven detonations are verified using a one-step viscous detonation model previously examined by Romick et al. (Reference Romick, Aslam and Powers2012) in the CJ limit. Here, a simple study was performed at several piston velocities, for a fixed non-dimensional activation energy of
$E_{a}=29.98$
which predicted a period-four detonation in the CJ case. As the focus of this examination is the energy composition driving the detonation at late times, only a domain near the front travelling at an average detonation velocity is considered. The domain is taken to be sufficiently large to encompass several half-reaction zone lengths. This domain has chemical energy entering from the unreacted ambient upstream conditions. Energy leaves the domain behind the front due mainly to thermal energy released from the reaction, as well as a small contribution from energy associated with viscous stresses. In addition, depending on the magnitude of the piston velocity, kinetic energy either propagates out of the domain or enters the domain emanating from the energy input from the piston. However, only the composition of energy within this domain is of interest in this case; thus, by integrating the individual amounts of these three (thermal, chemical and kinetic) energies as well as the sum over this domain at each time, the evolution of the energy composition driving the detonation at the different piston velocities can be obtained. The sum of the energies in the domain can be written as
$$\begin{eqnarray}E(t)=\int _{x_{1}}^{x_{1}+L_{d}}\left({\it\rho}c_{v}T+{\it\rho}\left(1-{\it\lambda}\right)q_{r}+{\it\rho}\frac{u^{2}}{2}\right)\text{d}x,\end{eqnarray}$$
where
$E(t)$
is the total energy per unit area in the domain,
$x_{1}$
the position of the left end of the domain,
$L_{d}$
the fixed length of the domain,
$c_{v}$
the specific heat at constant volume for the mixture in this one-step model,
${\it\lambda}$
the reaction progress variable and
$q_{r}$
the heat release per unit mass of the reaction. The total domain length chosen for integration is 50 steady half-reaction zone lengths with 35 steady half-reaction zone lengths behind the front at the peak detonation velocity in the late time behaviour. This choice allows for an almost complete reaction to occur within the domain at any period of the cycles predicted.
Figure 15 shows the fraction of the total energy each component carries for piston velocities of
$\overline{u}_{p}=0.0,0.5,1.0,1.1,1.2,1.3,1.4,1.5\times 10^{5}~\text{cm}~\text{s}^{-1}$
in the late time behaviour of the propagating detonation. In the ZND profile for this activation energy, the CJ particle velocity is
$0.96\times 10^{5}~\text{cm}~\text{s}^{-1}$
; therefore, it is anticipated piston velocities above this velocity would have an effect on the late time behaviour. The majority of the energy is carried in the thermal mode due to the energy released by the reaction, depicted in figure 15(a), with all piston velocities predicting a mean contribution near 75.5 % of the total energy. At sufficiently high piston velocities, a stable detonation is formed. This can be most easily seen in figure 15(b,c), where the chemical and kinetic modes are shown, respectively. In figure 15 the light grey dashed horizontal lines are for the highest piston velocities. For these stable detonations, the kinetic mode carries more energy than the chemical mode due the high piston velocities. However, as the final piston velocity is lowered, a period-one detonation propagates at late times as seen in figure 15 by the thick dark grey curve for
$\overline{u}_{p}=1.3\times 10^{5}~\text{cm}~\text{s}^{-1}$
and the solid light grey curve for
$\overline{u}_{p}=1.2\times 10^{5}~\text{cm}~\text{s}^{-1}$
. Furthermore, at the higher piston velocity
$\overline{u}_{p}=1.3\times 10^{5}~\text{cm}~\text{s}^{-1}$
, the kinetic and chemical modes are nearly identical in magnitude.
Fraction of total energy within the integration domain for (a) thermal, (b) chemical and (c) kinetic energies in the one-step model with
$E_{a}=29.98$
for several piston velocities.
At an even slower piston speed, a period-two detonation is predicted, depicted by the thinner dark grey curve. In this case,
$\overline{u}_{p}=1.1\times 10^{5}~\text{cm}~\text{s}^{-1}$
, the chemical mode is now dominant over the kinetic mode. Moreover, the ratio between the chemical and kinetic modes at its maximum is 2.24. The lowest three piston velocities examined,
$\overline{u}_{p}=0.0,0.5,1.0\times 10^{5}~\text{cm}~\text{s}^{-1}$
, exhibit period-four CJ behaviour. The separations between these CJ detonations are almost exclusively temporal offsets, and the curves would nearly be coincident if these offsets were removed. As the
$\overline{u}_{p}=1.0\times 10^{5}~\text{cm}~\text{s}^{-1}$
case is minutely above the CJ particle velocity for the ZND profile, there is only a weak contribution from the piston affecting the energy composition. These piston velocities are shown in figure 15 by the solid, dashed and dashed-dotted curves, respectively. For the CJ case, which is present for the
$\overline{u}_{p}=0.0,0.5\times 10^{5}~\text{cm}~\text{s}^{-1}$
cases, the maximum ratio between chemical and kinetic modes grows to 2.95 which is the maximum for this particular activation energy.
As discussed earlier, identifying the beginning and end of the limit cycle is challenging for multistep kinetics. However, given sufficiently long computation time it is presumed that fully relaxed limit cycles would be predicted. Therefore, it is hypothesized that these explanations which hold for one step kinetics also extend to multistep kinetics.
6. Conclusions
An investigation of one-dimensional piston-supported hydrogen–air detonations in the presence of mass, momentum, and energy diffusion has shown that as the final supporting piston velocity is decreased, the long time behaviour becomes more complex. For detonation propagating into a stoichiometric mixture at 1 atm and 293.15 K there are several distinct phases in the long time behaviour. At sufficiently high piston velocities, the detonation relaxes to a stable steadily propagating detonation wave. Below a critical supporting piston velocity, a pulsating detonation with a single high-frequency oscillation is predicted. In comparison with the inviscid analogue, the inclusion of physical diffusion delays this neutral stability boundary by
${\sim}2\,\%$
with respect to the supporting piston velocity. In addition, once the viscous analogue has begun pulsating, the addition of physical viscosity reduces the magnitude of oscillations by 40 % at times. This suggests that using methods in which numerical diffusion is important, such as filtering, large eddy simulation, or a calculation on a coarse grid, may be problematic in this regime as it may shift bifurcation points or reduce the amplitude of oscillations significantly. Lowering the supporting piston velocity further gives rise to a pulsating detonation in which both low and high frequencies exist and carry a significant amount of the pulsation energy at late times. Near this transition to the dual mode behaviour, the effect of adding viscosity has been weakened; with the magnitude of oscillations being reduced to less than 0.1 % difference between the inviscid and viscous cases. Below another critical point in supporting piston velocity, the long time behaviour becomes dominated by the low-frequency oscillations. In this regime, the inviscid and viscous cases are nearly indistinguishable due to the intrinsic instability growing in strength compared with the high-frequency mode. Within in this regime, a nearly period-doubling phenomenon is predicted giving rise to a period-two detonation. At sufficiently low piston velocities, the long time detonation pressure exhibits many different local maxima in time which is indicative of chaotic detonation.
In addition to examining the long time detonation pressure behaviour in the time domain, it was also examined in the frequency domain using the PSD to gain insight into how the energy variance changes versus supporting piston velocity. The fine resolutions used in this study enabled detailed results to be obtained from the harmonic analysis. It was found that in the high-frequency mode of pulsations the spectrum blue shifts towards lower frequencies as the supporting piston velocity was lowered with the fundamental frequency moving from
$3.41~\text{MHz}$
at
$\overline{u_{p}}=1.410\times 10^{5}~\text{cm}~\text{s}^{-1}$
to
$2.71~\text{MHz}$
at
$\overline{u_{p}}=1.320\times 10^{5}~\text{cm}~\text{s}^{-1}$
. After the second bifurcation point when the pulsations exhibit both a high- and a low-frequency oscillation, the spectrum begins to shift more energy to the lower frequency as the supporting piston velocity is lowered. In this dual mode region, there exists side banding around the high-frequency mode and its harmonics due to the interaction with the active low frequency mode. In addition, the spectrum continues to blue shift but at a slower rate than in the high-frequency mode. Though small, the addition of physical viscosity does have a stabilizing effect near the transition which is exhibited by the higher frequency carrying more energy in the viscous case than the inviscid analogue. As the supporting piston velocity is lowered further, the low-frequency mode becomes the dominant mode of the pulsations. In this region, the PSD spectrum eventually exhibits subharmonic frequencies indicating a near period-doubling phenomenon has occurred. At sufficiently low piston velocities, the long time behaviour has many active modes which is usually indicative of a chaotic detonation.
In addition, the form of initialization was briefly addressed and it was found that at long times the detonation relaxed to the same steadily travelling detonation. A limited validation study was conducted, and a frequency of
${\it\nu}=0.97~\text{MHz}$
was predicted for an overdriven detonation of
$f=1.10$
which is in agreement with the observation of 1.04 MHz by Lehr (Reference Lehr1972) in shock-induced combustion experiments around spherical projectiles in a hydrogen–air mixture at 0.421 atm and 293.15 K. Furthermore, the effect of the threshold parameter,
${\it\epsilon}$
, was examined, and it was found that only minuscule differences between the detonation pressure for a stable detonation indicating the WAMR method was in the convergent regime.
Acknowledgement
The authors thank the US Department of Energy for support of this work.
