A Zel’dovich–von Neumann–Döring-like detonation wave in a multi-temperature mixture

The detonation wave structure is analysed in a binary mixture undergoing a reversible chemical reaction represented by $A_{r}\rightleftharpoons A_{p}$ . It is assumed that the flow satisfies the proper basic assumptions of the Zel’dovich–von Neumann–Döring (ZND) detonation model, namely the flow is one-dimensional and the shock is represented by a jump discontinuity, but the assumption of local thermodynamic equilibrium is disregarded. This allows us to deeply investigate the coupling between the detonation structure of overdriven detonations and its chemical kinetics. The thermodynamic non-equilibrium effects are taken into account in the mathematical description, using the model of a multi-temperature mixture developed within extended thermodynamics, which has been proved to be consistent with a kinetic theory approach. The reaction rate is then enriched with terms that take into account the temperatures of the constituents. The results show that the temperature difference between components within the detonation wave structure, which describes thermodynamic non-equilibrium, is driven by the chemical reaction. Numerical computations confirm the existence of non-monotonic profiles in the reaction zone of overdriven detonations which are sensitive to changes in the activation energy and reaction heat.

Mathematical modelling of chemically reacting gas mixtures within continuum theories usually relies upon models that comprise the conservation laws of mass, momentum and energy of the mixture, and an additional balance law that describes the progress of the reaction by tracking the rate of change of product concentration, driven by the reaction rate. These models are, on the one hand, well founded in molecular gas theory. On the other hand, they are mathematically tractable and provide a good description of real processes. This classical approach is well explained in seminal works (Hirschfelder et al. 1954;Courant & Friedrichs 1977;Zel'dovich & Raizer 2002), as well as in more recent monographs devoted to combustion and detonation problems (Fickett 1985;Williams 1985;Lee 2008).
Versatile problems that appear in the detonation processes are both physically interesting and mathematically demanding. Physical aspects and various approaches that emerge in their study are well documented (Fickett 1985;Williams 1985;Lee 2008). On the other hand, the mathematical complexity of the detonation problem motivated the construction of model equations, like the one proposed by Majda (1981). It retained the basic features of the detonation process utilizing a reactive Burgers equation analogue. Such a model was analysed for the solution properties (Razani 2002), as well as for the stability of the detonations (Liu & Yu 1999;Szepessy 1999;Lyng & Zumbrun 2004). The existence of a solution for the complete set of detonation equations (Euler equations and reaction progress equation) was studied by Gardner (1983) and Hesaaraki & Razani (2001). A then up-to-date review, comprising the detonation structure, was given by Clavin (2004). Meanwhile, certain new aspects of detonation waves emerged within the continuum framework, e.g. radial propagation (Kasimov & Korneev 2014), multi-dimensional self-sustained detonations (Faria, Kasimov & Rosales 2015) and decay of plane detonation waves (Clavin & Denet 2018), so that detonation processes still remain in focus. Apart from the continuum approach, detonation problems were also analysed using the kinetic theory of gases as a starting point, and deriving macroscopic models as moment equations of distribution functions. In this approach the system was usually closed at Euler level (Conforto et al. 2003(Conforto et al. , 2004Carvalho & Soares 2012;Marques Junior et al. 2015), leading to hyperbolic systems of balance laws (the exception is Conforto, Monaco & Ricciardello (2014), where the closure was done at Navier-Stokes level). Finally, the intricacy of the numerical computation of detonation profiles due to the scale disparity brings additional difficulties. The computations should be carefully performed in order to obtain reliable results (see Bdzil & Stewart 2007;Cael et al. 2009).
The continuum models mentioned above can be enriched if one digs deeper into the structure of the medium. One of the simplest, and conceptually reasonable, ways is to take into account the information about the constituents by extending the list of state variables that describe the medium. In fact, this is the aim of the mixture theory developed within the framework of extended thermodynamics (Müller & Ruggeri 1998, pp. 79-104). In fact, in this approach, the state of the mixture is described by the fields of mass density, velocity and temperature of each constituent (Ruggeri & Simić 2007;Ruggeri 2009). This choice of state variables, and appropriate governing equations, enables more sophisticated modelling and acquiring more information about non-equilibrium processes. Moreover, this model for mixtures possesses two main features of extended thermodynamic models: it is dissipative, and it is hyperbolic. Nevertheless, it can be put within the classical framework either by means of Maxwellian iteration (Ruggeri & Simić 2009, or directly (Gouin & Ruggeri 2008).
1.1. Assumptions and the aim of the paper The aim of this paper is to propose a variant of the Zel'dovich-von Neumann-Döring (ZND) detonation model within the framework of multi-temperature extended thermodynamics of mixtures, and to investigate the detonation wave structure problem. In this approach the basic assumptions of the ZND model will be retained: the flow is assumed to be one-dimensional, the shock is described by a non-reactive jump discontinuity, and the chemical reaction occurs in a finite reaction zone following the shock front. Finally, the study will be concerned with overdriven detonations that propagate with a speed greater than the Chapman-Jouguet one. To that end, certain generalizations of the standard ZND model will be made. First, it will be assumed that the chemical reaction is reversible, restricting the analysis to the simplest possible mechanism, represented by A r A p , (1.1) with only two constituents -reactants (r) and products (p) of the forward reaction. This representation includes the symmetric reaction A 1 + A 1 A 2 + A 2 , which has been considered in previous works about detonation and linear stability (see Carvalho & Soares 2012;Marques Junior et al. 2015). Reaction (1.1) can also represent an inelastic transition of a polyatomic molecule from one energy state to another energy state and, in this case, A r and A p denote two different energetic components of the same gas. Second, the rather restrictive assumption of local thermodynamic equilibrium will be replaced with a set of new assumptions about the constituents and their dynamics, at both molecular and macroscopic scale: (i) the molecular masses are equal, m r = m p = m; (ii) the molecular diameters are equal, d r = d p = d; (iii) the binding (formation) energies of the molecules are different, r = p ; (iv) during elastic collisions, molecules have monatomic-like behaviour, and both constituents will be modelled as hard spheres; (v) the ratios of specific heats are equal, γ r = γ p = γ = 5/3; (vi) the macroscopic state of the constituents is described by the fields of mass density ρ α , velocity v α and temperature T α , with α = r, p; and (vii) the constituents are assumed to be Eulerian fluids, and thus the viscosity and heat conductivity of the constituents do not influence the detonation structure and are not taken into account.
The physical motivation for the choice of chemical reaction (1.1) and the set of assumptions (i)-(vii) is to capture, in the simplest possible setting, the scenario that is typical for overdriven detonations. Namely, in the analysis of pathological detonations, Sharpe (1999), and also Khokhlov (1989), mentioned that overdriven detonations have continuous but non-monotonic profiles, in which the turning point corresponds to the state with maximum heat release. So far, the simplest type of chemical reaction that allows this kind of behaviour is the two-step reaction A → B → C, where the first reaction is exothermic, while the second one is endothermic. In a recent paper, Marques Junior et al. (2015) used a reversible reaction of the kind (1.1), but did not comply with the expected behaviour of the state variables in overdriven detonations. Our aim is to bring into play the thermal non-equilibrium through the multi-temperature assumption, and to show that the standard overdriven detonation scenario is possible in (1.1) due to different temperatures of the constituents.
1.2. Discussion of the assumptions Assumptions (i) and (vi) are the most important ones for the generalizations in our model. Assumption (i), although it seems counter to our intentions for generalization, is a direct consequence of the mass conservation in chemical reaction (1.1). It induces a kind of singular behaviour of the mixture, since the molecules of the constituents cannot be distinguished in a mechanical sense due to their equal masses. This observation will be summarized as the basic assumption, and it will simplify the model -the difference of the velocities will be disregarded and a single-velocity model will be considered. On the other hand, assumption (vi) introduces the multi-velocity and multi-temperature structure into the model. This structure is already embedded into the rational thermodynamics of mixtures (see Truesdell 1984, pp. 219-236), and has become fundamental in the extended thermodynamic theory of mixtures (Ruggeri & Simić 2007) and the kinetic theory of reactive mixtures (Bisi, Martalò & Spiga 2011Murugesan, Sirmas & Radulescu 2017). Furthermore, it is of macroscopic kind, but at the same time it goes deeply into the structure of matter and resonates with molecular assumptions, thus reflecting the primary aim of extended thermodynamics -to bridge the gap between the macro-and meso-scales. In the context of chemically reactive mixtures, it aims to describe the internal energy transfer from reactants to products of the reaction, or vice versa, and the corresponding opposite kinetic energy transfer. Such an assumption has been proved to be reasonable when the constituents have large discrepancies in molecular masses, and thus considerably different internal energies (temperatures) (see Bird 1968;Kosuge, Aoki & Takata 2001;Madjarević, Ruggeri & Simić 2014). This observation came from studies in which non-reactive mixtures were analysed, and mutual interactions between the constituents were of purely mechanical nature. Our model will bring another novelty by applying the multi-temperature assumption to chemical interactions as well: the reaction rate for the reversible chemical reaction will be enriched with terms that take into account the constituents' temperatures. It will be shown that this generalization, in conjunction with large initial discrepancy of concentrations of the constituents, will be sufficient to drive the constituents' temperatures apart within the detonation wave. In such a way, we shall be in a position to prove that the difference of the temperatures, driven by the chemical reaction, influences the structure of the detonation profile.
Although our study is concerned with the relation of the multi-temperature assumption to the detonation profile, it does not have the aim to cover all the involved aspects and possible consequences. It will be shown in the discussion that it can be legitimately applied to multi-step chemical reactions or to other complex chemical processes, and used to describe both the structure of pathological detonation waves (Fickett & Davis 1979;Sharpe 1999) and the effects resulting from the competing exothermic and endothermic reactions. In particular, in conjunction with reversible chemical reactions, it predicts the appearance of the non-monotonic profiles expected in the case of overdriven detonations, which we shall deal with in this paper.
Furthermore, the multi-temperature assumption may be related not only to different kinetic temperatures of the constituents, but also to different degrees of freedom of polyatomic molecules through vibrational temperature, as well as to different phases.
In particular, with regard to the multi-temperature formalisms applied to treat detonation, we quote, in particular, the relevant work about the modelling, structure and other aspects of detonation solutions in two-phase reactive flows (Powers, Stewart & Krier 1990a,b) and the important contributions about the non-equilibrium effects in gaseous detonation systems resulting from the interplay between vibrational relaxation and chemical reactions (Taylor, Kessler & Oran 2013;Shi et al. 2017).
1.3. Structure of the paper The paper is organized as follows. In § § 2 and 3 we shall provide a proper background for the modelling of ZND detonation waves. Section 2 will review the basic aspects of the multi-temperature model for a binary reactive mixture, while § 3 will explain the main features of the steady ZND detonation waves in the context of reversible chemical reactions. Section 4 will contain the formulation of the steady detonation problem and the reduction of the model based upon the basic assumption regarding diffusion. After that, qualitative aspects of the reacting flow profiles will be analysed, which include the equilibrium state and the existence of non-monotonic profiles. In § 5 continuous multi-temperature steady reacting flow profiles are computed numerically for overdriven detonations, and compared with appropriate single-temperature profiles. Since the solution of the detonation problem depends on several parameters, their influence will be thoroughly analysed. It will be shown that, for certain values of the parameters, the multi-temperature model yields qualitatively different detonation profile in comparison to the single-temperature one. For other values of the parameters, the qualitative picture is the same. Nevertheless, the influence of two temperatures in the model is inevitable in the structure of the profiles. The paper ends with a summary and some possible studies to be considered as a continuation of our paper.

The multi-temperature model of reactive binary mixture
The framework for the detonation wave structure problem with reversible chemical reaction in a binary mixture will be the multi-temperature (MT) model developed within extended thermodynamics (Ruggeri & Simić 2007). A similar model was recently applied to the description of flame structures (Brini 2009). The main feature of the model is the introduction of the velocities and the temperatures of the constituents, with the aim of capturing non-equilibrium effects in the simplest possible constitutive framework. Although it may seem quite unusual to extend the list of state variables of a macroscopic model in such a way, it resonates with the appropriate models derived within the kinetic theory of gases (Goldman & Sirovich 1967), and thus goes along with the main intention of rational extended thermodynamics (Müller & Ruggeri 1998) to bridge the gap between the macroscopic and mesoscopic levels of description. The rationale for this approach in detonation problems comes from the fact that chemical reactions are processes that could drive the medium out of equilibrium, in a chemical sense, as well as in mechanical and thermal senses.

Governing equations
The main modelling assumption is that each constituent of the mixture obeys the same balance laws as a single fluid, but mutual interaction of the constituents must be taken into account (Truesdell 1984). The governing equations for an n-component mixture (α = 1, . . . , n) then read (see Ruggeri & Simić 2007): where ρ α , v α and ε α are the mass density, the velocity and the internal energy of the α-constituent, T α and q α are the stress tensor and the heat flux, while the source terms τ α , m α and e α take into account the mutual interactions of the constituents. Moreover, v α = |v α | and similar notations are adopted for other velocities from here on. It is an axiom of the mixture theory that conservation laws of mass, momentum and energy of the whole mixture must hold, imposing the following constraints on the source terms: In accordance with our assumption (vii), the stress tensor and the heat flux will be reduced to T α = −p α I, q α = 0, (2.3a,b) where I is the unit tensor. For the sake of brevity we shall omit the general analysis of the MT model of mixtures, which can be found elsewhere (Ruggeri & Simić 2007;Madjarević & Simić 2013;Madjarević et al. 2014). Instead, the case of a binary reactive mixture will be presented. The governing equations consist of conservation laws of mass, momentum and energy for the mixture, and balance laws of mass, momentum and energy for one constituent (Ruggeri & Simić 2007). Since it is common practice in simple modelling of chemical reactions to monitor the behaviour of the product of the forward reaction, we decided to include the balance laws for the product in the governing equations, so that they read: and ∂ ∂t (2.9) Here the following quantities are introduced: ρ = ρ r + ρ p mass density of the mixture, z = ρ p /ρ concentration of the products, v = 1 ρ (ρ r v r + ρ p v p ) mean velocity of the mixture, p = p r + p p total pressure of the mixture, J = ρ p u p = −ρ r u r diffusion flux vector, u α = v α − v (α = r, p) diffusion velocities, ε I = 1 ρ (ρ r ε r + ρ p ε p ) intrinsic (thermal) internal energy density of the mixture, ε = ε I + 1 2ρ (ρ r u 2 r + ρ p u 2 p ) internal energy density of the mixture, where g α are dynamic enthalpies. To close the system of equations (2.4)-(2.9), one has to introduce constitutive assumptions about pressure and internal energy density, as well as to determine the structure of the source terms. The latter is dictated by the general principles of extended thermodynamics -Galilean invariance and the entropy principle. Galilean invariance determines the velocity dependence of source terms (Ruggeri & Simić 2007): whereτ p ,m p andê p denote velocity-independent mass, momentum and energy production rate, respectively. Information about the chemical reaction is contained in the source term τ p of the mass balance law (2.7). Nevertheless, it is reflected in all other source terms as well, due to Galilean invariance (2.12). The entropy principle imposes a further restriction on the source terms by requiring that every thermodynamic process is compatible with an entropy inequality. Its consequences will be discussed in the sequel.
2.2. Constitutive assumptions and the state variables Since viscosity and heat conductivity are not taken into account in our model, the constitutive assumptions are related to thermal and caloric equations of state. We shall restrict our study to a binary mixture of ideal gases that consists of reactants and products. As we are analysing the simplest possible model of a reversible chemical reaction, A r A p , according to assumption (i) the masses of the constituents ought to be equal due to mass conservation, m r = m p = m. Consequently, the equations of state have the following form: where k is the Boltzmann constant, c vα (T α ) are the specific heats at constant volume and α are the binding energies. Assuming c vp = c vr = c v = const., and r = c v T 0 , the internal energies take the form: . The temperatures of the constituents are not equal in general, and we have to introduce the average (mean) temperature of the mixture as a global variable. Following Gouin & Ruggeri (2008) and Ruggeri & Simić (2009), we shall define the average temperature using the definition of the intrinsic part of the internal energy density: where we have taken into account the relations ρ p = ρz and ρ r = ρ(1 − z). It is also convenient to introduce the difference of temperatures, as a state variable, which also serves as a measure of thermal non-equilibrium. In such a way, the temperatures of the constituents can be expressed in terms of z, T and Θ, as follows: The total pressure p and the intrinsic internal energy density ε I of the mixture can be expressed in terms of T: where we have taken into account assumption (v), that is, γ p = γ r = γ . Having in mind the structure of the governing equations (2.4)-(2.9), we may fix the set of state variables. Although the first choice is usually (ρ, v, T, ρ p , v p , T p ), using the relations given above, we shall adopt U = (ρ, v, T, z, J, Θ).

The source terms
In our model, the source terms describe the mutual interactions between the constituents. Their structure, compatible with entropy inequality, was determined by Ruggeri & Simić (2007) when a non-reactive mixture is considered. However, the extended thermodynamics has its limits and the phenomenological coefficients cannot be determined within its framework. Instead, one has either to take into account experimental results, or to make a link with a more refined approach (e.g. kinetic theory) to determine them (see Madjarević & Simić 2013;Madjarević et al. 2014). We shall proceed in the latter way.

The reaction rate
For the reactive mixture considered in this paper, the reaction rate is obtained from the kinetic theory; more precisely, from the Boltzmann equation for chemically reactive mixtures, assuming reactive differential cross-sections with activation energy, and a particular form of the velocity distribution functions in order to obtain an approximate solution to the reactive kinetic equations. Since we are considering an MT mixture, an appropriate form for such an input function should be a Maxwellian distribution in the constituent rest frame, where c α = |c α |, c α being the molecular velocity of each constituent, and n α = ρ α /m the constituent number density. However, the evaluation of the collision integrals applying to the input functions (2.19) becomes a very difficult task and the computations lead to intractable expressions. To overcome such difficulty we approximate function (2.19) through an expansion at the first order around a Maxwellian in the mixture rest frame with temperature T and mean velocity v: Here (v αi − v i ) and (c αi − v i ) represent the spatial components of the diffusion velocity and of the peculiar velocity, respectively, and ∆ α = (T α − T)/T is the normalized temperature difference. The input function (2.20) represents a linearization of the Maxwellian distribution (2.19), with respect to the diffusion velocity and temperature difference. Adopting reactive differential cross-sections with activation energy, and taking into account the indistinguishability of the colliding particles, the reaction rate for the mass concentration z of the products results in (Kremer, Bianchi & Soares 2007) where f and b denote the forward and the backward activation energies of the chemical reaction. The reaction rate (2.21) shows a second-order form, in the sense that it depends explicitly on the square of the product concentration z and of the reactant concentration (1 − z). This is a direct consequence of the fact that it was derived from the Boltzmann equation with binary reactive encounters. A detailed account of other aspects of the mathematical modelling of reactive flows may be found in Giovangigli (1999). The thermodynamic equilibrium condition, chemical and thermal, corresponds tô τ p = 0 and, therefore, we will have thermodynamic equilibrium if is the reaction heat. The last condition in (2.22) represents the law of mass action of the model and gives a constraint on the chemical equilibrium values of the product concentration z and mixture temperature T.

The momentum production rate
The momentum production rate in a binary mixture has the form: where ψ pp is a phenomenological coefficient related to the relaxation time τ D for diffusion. By linearizing the source term in the neighbourhood of the local equilibrium state, and taking into account the definition of the diffusion flux, one obtainŝ (2.24) 2.3.3. The energy production rate The energy production rate in a binary mixture readŝ where θ pp is a phenomenological coefficient related to the relaxation time τ T for the temperature. By linearizing the source term in the neighbourhood of the local equilibrium state, and using the difference of temperatures Θ as non-equilibrium variable, one obtainsê In expressions (2.24) and (2.26), the relaxation times τ D and τ T , which are part of the phenomenological coefficients, cannot be determined by the methods of extended thermodynamics. It was shown in Madjarević & Simić (2013) and Madjarević et al. (2014), using the arguments of kinetic theory, that the relaxation times for diffusion and temperature can be related to diffusivity of the binary mixture. In our notation, the diffusivity D for the hard-sphere model reads: ( 2.27) Taking into account that m p = m r = m, the relaxation times reduce to the following simple form: Note that τ T > τ D , which is a typical result of kinetic theory.
2.4. The basic assumption and reduced form of the governing equations Equations (2.4)-(2.9), with the source terms given by (2.21), (2.24) and (2.26), constitute a complete set of governing equations for a chemically reacting MT binary mixture. In their derivation we exploited all the simplifying assumptions mentioned in the introduction. However, it was mentioned that assumption (i) implies further simplifications. It is appropriate at this stage to introduce the basic assumption that will imply the reduction of the system of governing equations, and simplification of the remaining ones.

The basic assumption
Since reactants and products have equal molecular masses, m r = m p = m, the molecules of the constituents are indistinguishable in a mechanical sense. Therefore, the diffusion of the constituents can be ignored and the momentum balance law (2.8) can be dropped from the list of governing equations.
The basic assumption is a consequence of the simplicity of the chemical reaction (1.1). The fact that the masses of the constituents are equal implies that one cannot distinguish the molecules of reactants and products during mechanical collisions. As a consequence, the mean (macroscopic) velocities of the constituents are equal, v r = v p = v, and the diffusion velocities vanish, u r = u p = 0, as well as the diffusion flux, J = 0. However, the mean energies of the constituents may not be equal, which implies thermal non-equilibrium, T r = T p , and chemical reaction may play a role in their approach to equilibrium. The MT assumption is in this case motivated by the presence of chemical reaction, rather than the difference of masses as in previous studies (

The reduced form of the governing equations
As a consequence of the basic assumption, the set of governing equations reduces to five equations, namely (2.4)-(2.7) and (2.9), for the state variables U = (ρ, v, T, z, Θ), with production termsτ p andê p given by expressions (2.21) and (2.26).
The reduced set of equations describes the evolution of the considered binary reactive mixture, in a regime of thermal and chemical non-equilibrium in which diffusion, heat conduction and viscosity are not taken into account. Such a set of equations is the one considered in the sequel to define the mathematical setting for the determination of the detonation wave solutions.

Steady ZND detonation waves
In this section we follow the classical ZND model, advanced by Zel'dovich, von Neumann and Döring (see Fickett 1985;Williams 1985;Lee 2008), for a one-dimensional detonation wave, and propose its generalization taking into account the reversibility of the chemical reaction and the MT assumption.
The configuration of the ZND detonation wave consists of a plane, non-reactive shock wave propagating with constant speed s in the direction orthogonal to the singular surface, followed by a finite reaction zone where the chemical reaction takes place and eventually reaches equilibrium.
The spatial structure of the steady detonation wave is determined by solving two different problems in a convenient parametric space. The first problem is purely algebraic and consists in determining the post-shock state -the von Neumann state N -as the solution to the non-reactive jump Rankine-Hugoniot conditions. The second problem is the initial value problem for the determination of the continuous reacting flow in the reaction zone, with initial data at the von Neumann state. The final state S, at the end of the reaction zone, can be obtained either from the second problem, as the state of the continuous reacting flow for which the reactive mixture reaches the chemical and thermal equilibrium, or from a Rankine-Hugoniot analysis combined with the equilibrium condition, chemical and thermal.
In what follows we refer the governing equations to the frame attached to the shock wave, and establish the analytic part of the modelling.

Governing equations in the moving frame
According to the ZND model, we shall assume that the wave moves in the x-direction, from left to right, and that the flow is one-dimensional, v = (v, 0, 0). The wave configuration is supposed to be steady in the shock-attached frame, and we introduce the steady variable ξ and the relative velocity u of the mixture with respect to the shock front: In the steady frame attached to the shock wave, the reaction zone extends from x 0 to x S , where x 0 represents the location of the shock front. The state just behind the shock, located at x = x 0 , is the von Neumann state N, where the chemical reaction is triggered, and the one located at x = x S , at the end of the reaction zone, is the final steady state S, where the mixture reaches chemical and thermal equilibrium. Ahead of the shock front, that is, for x > x 0 , the quiescent mixture is at rest in its initial state I, where the rate of the chemical reaction is completely negligible. Assuming steady, one-dimensional travelling wave structure of the detonation wave, in which the vector of state variables U = (ρ, u, T, z, Θ) depends on a single independent variable, say U(ξ ) = U(x − st), the governing equations are transformed into a set of ordinary differential equations: The source termsτ p andê p are used in linearized form (2.21) and (2.26), respectively. Formally, the governing equations of the standard ZND detonation waves are (3.2)-(3.5). Here, they are extended with (3.6), which describes the profile of the temperature difference Θ.

3.2.
Algebraic analysis of the detonation wave As the first step, the algebraic part of the problem will be studied. It will outline the basic structure of the detonation wave, including the important notion of equilibrium Hugoniot curve, peculiar for reversible chemical reactions.

Rankine-Hugoniot conditions
The governing equations of conservative type, equations (3.2)-(3.4), lead to the well-known algebraic Rankine-Hugoniot (RH) conditions that relate a state behind the shock wave to the initial state ahead the shock. Such conditions connect the fluxes of the conserved quantities across the wave front, and are given by where ρ 0 and p 0 are the mass density and pressure at the unperturbed initial state. Moreover, z 0 = 0 and Θ 0 = 0 at the initial state. Additionally we have assumed that the initial state is at rest, so that v 0 = 0 and u 0 = −s. For each fixed value of the detonation velocity s, the corresponding von Neumann and equilibrium states can be determined from the RH conditions (3.7)-(3.9), with additional proper conditions, as will be explained in the sequel. Introducing the specific volume of the mixture, ν = 1/ρ, the states (ν, p) in the reaction zone, compatible with a strong detonation wave (Lee 2008), are characterized as a function of the product concentration z, with the shock speed s as parameter, by where Υ (z) is defined by (3.11) In the derivation of (3.10) we used the state equations (2.18); see Marques Junior et al. (2015) for a detailed analysis.

Von Neumann state
The von Neumann state is the state just behind the shock wave. Since the passage of the mixture through the shock front is a non-reactive process, the progress variable remains constant through the shock and the temperature difference does not change. Accordingly, z = 0 and Θ = 0 at the von Neumann state. Therefore, the state N can be obtained as the non-trivial solution to the RH conditions (3.7)-(3.9) complemented by the conditions z = 0 and Θ = 0. Equivalently, the state N can be obtained directly from conditions (3.10), together with the vanishing conditions for z and Θ, that is

Continuous reacting flow
The reacting flow in the reaction zone describes the continuous evolution of the mixture from the von Neumann state N, say at ξ = 0, to the final state S, say at ξ = ξ S , due to the chemical reaction and thermal equilibration. From the mathematical point of view, the ordinary differential equations (3.2)-(3.6) have to be solved for the fields ρ, u, T, z and Θ, with initial conditions at the state N.

Equilibrium final state
The final state S at the end of the reaction zone represents the state for which the reactive mixture reaches chemical and thermal equilibrium. It can be determined from expressions (3.10) with z and Θ satisfying the equilibrium conditions (2.22). Therefore, S is characterized by expressions (3.10) together with where we have used Υ (z) defined by (3.11). Observe that the second condition in (3.13) represents the law of mass action given in (2.22) parametrized by the shock velocity s. Additionally, since the state S represents the end point of the continuous reacting flow, in which the mixture reaches equilibrium and all state variables become constant, it can be obtained by integrating equations (3.2)-(3.6) from ξ = 0 to ξ = ξ S for which all the production rates vanish, as predicted by the equilibrium conditions (2.22) or (3.13).

Chapman-Jouguet state
The Chapman-Jouguet state (CJ) is a particular final state obtained for the minimum admissible value of the shock velocity s, say s j . Since it is an equilibrium final state, equilibrium conditions (3.13) must be satisfied. Therefore, the CJ state is the final state of a strong detonation for which p/p 0 is minimum and ν/ν 0 is maximum, see (3.10). This corresponds to having Υ (z) = 0 (see Marques Junior et al. 2015). Accordingly, we determine the CJ velocity s j and the corresponding product concentration z j by solving the equilibrium conditions (3.13), together with condition Υ (z) = 0. The CJ state is then obtained from (3.10) for z = z j and s = s j .

Hugoniot diagram
The detonation solution can be represented in the (ν, p)-plane through the Hugoniot diagram (Fickett 1985;Lee 2008;Conforto et al. 2014;Marques Junior et al. 2015) straightly connected to the RH conditions (3.7)-(3.9). As shown in figure 1, this is a qualitative representation of the solution in terms of Rayleigh lines and Hugoniot curves for fixed product concentration, as well as in terms of the equilibrium Hugoniot curve, which constitutes an important ingredient when the chemical reaction is reversible and the reactive mixture reaches a final equilibrium state (Lee 2008;Marques Junior et al. 2015).
For fixed values of the material properties, say the specific heat ratio γ , the binding energy difference and the activation energy f , each Rayleigh line gathers all states obtained for a fixed value of the wave velocity s corresponding to mass and momentum conservation. Each Hugoniot curve for fixed product concentration collects all states along which the total energy of the mixture is preserved. On the other hand, the equilibrium Hugoniot curve accommodates all final equilibrium states for fixed values of the material properties. The necessity for the equilibrium Hugoniot curve comes from the fact that we analyse a reversible chemical reaction. Such a curve is the locus of all thermal and chemical equilibrium states, parametrized by the wave velocity s.
In the Hugoniot diagram, we highlight two items, namely the CJ state, corresponding to the tangency point of the Rayleigh line to the equilibrium Hugoniot curve, and the equilibrium final state (ν eq , p eq ) for a fixed value of s (s s j ), obtained as the intersection point of the corresponding Rayleigh line with the equilibrium Hugoniot curve.
For the considered value s > s j , we also plot in figure 1 the fixed-concentration Hugoniot curve obtained for the corresponding value z = z eq (solid grey curve). The equilibrium state (ν eq , p eq ) for such s > s j , point S in figure 1, is then the intersection of three curves: the corresponding Rayleigh line, the equilibrium Hugoniot curve, and the z = z eq Hugoniot curve. Although in figure 1 it seems that the latter two curves almost coincide, they do not and they must be distinguished.

Continuous reacting flow -qualitative analysis
Certain features of the MT ZND detonation waves can be revealed by a qualitative analysis. On the one hand, a hierarchy of models can be established and the MT model gradually simplifies to the standard ZND model. On the other hand, certain properties of the solution, peculiar for MT models, can also be recognized.

Steady detonation equations and hierarchy of the models
The governing equations (3.2)-(3.6) will be written in dimensionless form using the following scaled variables (subscript 0 refers to the unperturbed initial state in front of the shock wave): (4.1a−g) where the reference time τ 0 is given by Although the choice (4.2) is taken for convenience, τ 0 and relaxation times τ D and τ T (2.28) can be related to the so-called collision time scale, i.e. the mean free time of particles between collisions in equilibrium (Bardos, Golse & Levermore 1993). Omitting, for brevity, the tilde in the quantities (4.1), the dimensionless governing equations of the steady detonation are as follows: The governing equations consist of the conservation laws (4.3)-(4.5) of mass, momentum and energy for the mixture, and the balance laws (4.6) and (4.7) of mass and energy for the product. The balance law of momentum for the product is dropped from the set of governing equations in accordance with the basic assumption.

Governing equations solved with respect to first derivatives
The governing equations (4.3)-(4.7) can be solved with respect to first derivatives: (4.14) This form will be useful for the analysis of certain qualitative aspects of the solutions.

Hierarchy of the models
In § 5 we shall provide a comparative analysis of the detonation profiles obtained by the solutions to (4.3)-(4.7), with those obtained through the solutions to the simpler models. This analysis will be achieved through the hierarchy of models obtained through the simplification of the reaction rate (4.8). Namely,τ p can be split into a single-temperature (ST) partτ ST p , which does not depend on Θ, and an MT partτ MT p , which depends on Θ, as follows: In such a way the following hierarchy of the models can be obtained and analysed in the sequel: (MT) the MT model is described by the system (4.3)-(4.7), with the complete reaction rate term (4.8); (qMT) the quasi-MT model consists of the same equations (4.3)-(4.7) as for the MT model, but with simplified reaction rateτ p =τ ST p that does not take into account the contribution of the different temperatures; and (ST) the ST model consists of (4.3)-(4.6), assumes that the reaction rateτ p =τ ST p and assumes a common temperature for both constituents throughout the whole process -this model, along with reversible chemical reaction, corresponds to the model studied in (Marques Junior et al. 2015).

4.2.
Properties of the solution The analysis of MT detonation waves possesses certain particular features which distinguish it from the standard analysis based upon reactive Euler equations. Such features will be qualitatively analysed in the sequel.

Equilibrium state
As mentioned above, the MT model of the steady detonation wave describes a continuous solution that connects the von Neumann state U N with the equilibrium one U eq . Thus, we analyse the boundary value problem for (4.3)-(4.7) in a (formally) unbounded domain ξ ∈ [0, ∞), with the following boundary data: 19a,b) where ξ S denotes the value of the steady variable at which the process terminates. Note that condition (4.19) leads to a necessary condition for equilibrium.
STATEMENT 1. If the system (4.3)-(4.7) is in an equilibrium state characterized by relations (4.19), then the state variables obey the following relations: 1 − z eq z eq = exp − T eq and Θ eq = 0. (4.20a,b) Proof. The terminal condition (4.19) implies that production terms in governing equations (4.6) and (4.7) vanish in equilibrium. In particular, the equilibrium condition for (4.6) implies that the reaction rate vanishes,τ p = 0, which provides a relation between equilibrium values z eq , T eq and Θ eq . However, the equilibrium condition for (4.7) requires simultaneous vanishing of the reaction rateτ p and the temperature difference Θ, which leads to Θ eq = 0. As a consequence, we have Θ eq = 0 ⇒τ MT p (z eq , T eq , Θ eq ) = 0, (4.21) so that the ST partτ ST p of the reaction rate vanishes independently, implying the law of mass actionτ ST p (z eq , T eq , ) = 0 ⇒ 1 − z eq z eq = exp − T eq , (4.22) which completes the proof.
Having in mind the analysis developed in § 3, it becomes obvious that the terminal state of the system has to be the equilibrium one, i.e. U(ξ S ) = U S = U eq , where the equilibrium state is uniquely determined by the intersection of the Rayleigh line and the equilibrium Hugoniot curve, together with the independent condition Θ S = Θ eq = 0. The additional condition (4.20b) is peculiar for the MT model and has to be considered with particular attention, as will be shown in the sequel.

Non-monotonic detonation profiles
Further qualitative analysis will be based upon twofold characterization of the equilibrium state, given by the terminal condition (4.19). A distinguished role in this analysis will be played by the temperature difference Θ. To facilitate the formulation of the forthcoming statements, we shall introduce the following notation for the state variables: In the case of reversible chemical reaction, typical ZND detonation profiles are monotonic (Marques Junior et al. 2015), except in the neighbourhood of the Chapman-Jouguet solution, the solution that connects the von Neumann state N with the Chapman-Jouguet state CJ. However, equilibrium conditions (4.20) of the MT model imply the possibility of occurrence of non-monotonic profiles. The existence of such profiles will be demonstrated first through a qualitative analysis. It will then be confirmed by numerical computation in the next section.
Proof. If z C and T C satisfy the equilibrium condition (4.20a), i.e. the law of mass action, then ρ C and u C are uniquely determined by (3.10) and (3.7) for the fixed value of s. Moreover,τ ST p (z C , T C , Θ C ) = 0. In case (i) it is obvious that the equilibrium state is reached, since both equilibrium conditions (4.20) are satisfied. In case (ii), Θ(ξ C ) = 0 impliesτ p (ξ C ) =τ MT p (z C , T C , Θ C ) = 0. Consequently, by a simple inspection of the system (4.10)-(4.14), one may conclude that dU(ξ C )/dξ = 0, which completes the proof.
The crossing condition has a simple interpretation. If the state variables U ST satisfy the law of mass action (4.20a), and lie on the Rayleigh line and the equilibrium Hugoniot curve simultaneously, the ST part of the reaction rate will vanish,τ ST p (ξ C )=0. However, if Θ(ξ C ) = 0, the process is not terminated, since the derivatives of the state variables do not yet vanish, and the curve U ST (ξ ) only crosses the equilibrium values at ξ C . This motivates the denomination crossing condition. Its importance lies in the fact that it may occur for certain values of the parameters, in contrast to the standard ST case.
As a conclusion of the analysis of non-monotonic profiles, one fundamental remark has to be given. Either in the case of crossing the equilibrium state (statement 2), or in the case of local extremum (statement 3), the main driving agent for continuation of the process is the temperature difference of the constituents.

Quasi-multi-temperature behaviour
In § 3 we defined the qMT model as the one consisting of (4.3)-(4.7), or equivalently (4.10)-(4.14), but with a simplified reaction rateτ p =τ ST p . Note that such a reaction rate does not depend on Θ, but only on the U ST state vector. This fact has a simple consequence that will be exploited in our numerical computations.
STATEMENT 4. Assume that the reaction rate has the formτ p =τ ST p . Then (4.3)-(4.6), or equivalently (4.10)-(4.13), have the same structure as in the ST case and completely determine the solution for U ST (ξ ). Equation (4.7), or equivalently (4.14), is decoupled from the rest of the system and serves for the determination of Θ(ξ ) once the ST part of the model is solved.
Proof. It is sufficient to note that the differential part of the system (4.3)-(4.6) has the same structure as in the ST case (Marques Junior et al. 2015), and thus does not depend on Θ. If moreoverτ p =τ ST p , then this part of the system has the ST structure, independent of Θ, the same as the one in Marques Junior et al. (2015). Therefore, the complete system (4.3)-(4.7) decouples into an ST part (4.3)-(4.6), which determines the behaviour of U ST (ξ ), and an MT part (4.7), which tracks the evolution of Θ(ξ ). The proof is then complete.

Continuous reacting flow profiles
The profiles of a continuous reacting flow can be determined numerically, by solving an appropriate initial value problem, as discussed in § 3.2. Namely, the governing equations of the problem are (4.3)-(4.7), for ξ ∈ (0, +∞), and the initial data are taken as where U N stands for the von Neumann state. By M 0 = u 0 /c 0 we denote the Mach number in the unperturbed initial state in front of the shock, with c 0 = √ γ p 0 /ρ 0 being the adiabatic speed of sound. It will be convenient to express the Mach number in an alternative form, for the sake of comparison with the results presented in Marques Junior et al. (2015): where s j is the dimensionless steady detonation wave velocity at the CJ state, and f is the so-called overdrive degree. For reasons that will become apparent in the sequel, a subsequent analysis of the reacting flow profiles will be restricted to overdriven detonations, f > 1, for which s > s j and the equilibrium state lies on the strong branch of the equilibrium Hugoniot curve (see figure 1). To avoid singularities in numerical computations, we have to adopt a small value of the initial concentration of the product at the von Neumann state. Thus, we shall assume z(0) = 0.01 in all the computations. Other possible singularities do not interfere with numerical simulations. In particular, we have ρ > 0 and u > 0, and thus ρ and u do not create any singularity; moreover, z < 1, since the chemical reaction is reversible, and no other singularity arises from z. Finally, since we analyse only the overdriven detonations, the flow in the reaction zone is subsonic and u 2 < γ T throughout the flow.
The latter singularity, corresponding to the sonic condition u 2 = γ T, is rather delicate and has to be treated with care. When the detonation process comprises an exothermic reaction followed by an endothermic one, there is a possibility for a pathological detonation to occur (Fickett & Davis 1979). There is a point within the profile where the flow becomes sonic, and the detonation is unsupported and ends up on the weak branch of the equilibrium Hugoniot curve (see § 3.1 in Sharpe (1999) for a comprehensive analysis). In this case, the flow described by the model (4.3)-(4.7) will reach the singularity at the sonic point. For that reason, the analysis is restricted to overdriven detonations, while the other possibilities will be outlined in the final section.
The terminal, i.e. the equilibrium, state is determined by condition (4.19) for ξ → ξ S . In the numerical computations, we perform the analysis on a finite domain ξ ∈ [0, ξ * ], ξ * ξ S , taking ξ * in such a way that the state variables reach the neighbourhood of the equilibrium state with sufficient accuracy. We shall not go for the determination of the position ξ S of the equilibrium state, since it can be determined only approximately and the result depends on the accuracy proposed beforehand.
The basic modelling novelty, i.e. the introduction of multiple temperatures, influences neither the von Neumann state nor the terminal equilibrium state. It may be reflected only in the reacting flow profiles. This is a consequence of the fact that von Neumann and equilibrium states are determined through the solution of the appropriate RH conditions and the vanishing of production terms, respectively, which both include the vanishing of the temperature difference, as described in § 3.2. The aim of the numerical computation of the reacting profiles is to analyse the influence of the temperature difference on its structure. It is important to note that the model (4.3)-(4.7) together with initial data (5.1) contains three parameters: the forward activation energy f , the overdrive degree f , and the binding energy difference (or the reaction heat Q = −2 ). Our aim is to analyse the influence of these parameters on the structure of the continuous reacting flow, and to compare these results with those presented in Marques Junior et al. (2015), i.e. with ST profiles. Note that in view of statement 4, there is no need for a comparison with qMT profiles, since they are the same as ST ones. In figures 2-6, the steady variable ξ is given on a logarithmic scale and the state variables u, p and T are given in units of their corresponding values u N , p N and T N at the von Neumann state. 5.1. Case 1: Influence of f on the profile First we consider the influence of the activation energy of the forward reaction f on the structure of the continuous reacting flow profile. This is the case in which the MT assumption strongly reveals its influence and leads to completely new results in a qualitative sense. It is well known (Lee 2008) that the activation energy is strictly related to the temperature sensitivity of the chemical reaction. To compare our results with (Marques Junior et al. 2015), we adopted the same values of the parameters: Numerically computed profiles are presented in figure 2, for both MT and ST models.
For lower values of f , one can observe two qualitatively different types of profiles. The first one ( f = 2.0) is non-monotonic -it contains a crossing point with the state U ST eq , as indicated by statement 2, as well as at least one local extremum, as described in statement 3. This does not occur in the ST case at all. Moreover, it can be observed that the MT profiles are driven out of the von Neumann state faster than the ST ones, and this fact can be attributed to the presence of MT partτ MT p in the reaction rate.
When the activation energy is increased ( f = 10.0), the profiles become monotonic, like in the ST case. There is still, however, a difference between the MT and ST profiles: the MT profiles are driven out of the von Neumann state faster than the ST ones, but the ST profiles are steeper than the MT ones. Obviously, there exists a threshold value of the activation energy f which separates the non-monotonic profiles from monotonic ones. This issue will be tackled in the sequel.
As in Marques Junior et al. (2015), for very large values of the activation energy ( f = 30.0), the reaction rate is considerably smaller (see (4.8)) and the profile is much wider -the equilibrium state is reached at ξ eq ≈ 10 6 . Both MT and ST profiles are monotonic and have a similar behaviour as in the previous case. However, the difference between the state variables in the MT and ST cases is much less pronounced, e.g. z MT (ξ ) − z ST (ξ ) < 10 −3 throughout the profile. For these reasons, the profiles obtained for f = 30.0 are not displayed in the whole range of ξ , but just in the range that is significant for the first two cases. The analysis of the temperature profiles of the constituents sheds new light on the MT phenomenon. In the case of non-reactive mixtures, it is usually related to the fact that constituents have disparate masses. In our problem the masses of reactants and products are equal by assumption, m r = m p = m. Nevertheless, there is a temperature difference that cannot be ignored. A careful examination of the profiles, presented in figure 3, reveals that there appears an immediate temperature drop of the product temperature right behind the von Neumann state. This result suggests that the main reason for sustained thermal non-equilibrium between the constituents is not the mass difference, but the chemical reaction and the huge initial difference in the concentration of the constituents.
As a final remark, we would like to compare the non-monotonic profiles that appear in our model with the non-monotonic profiles in the standard model for ZND detonation waves. Namely, in the standard model, the non-monotonic temperature profiles appear for large values of the activation energy, and they are preceded by a longer induction zone and have steeper profiles (see Lee 2008, pp. 77-78). Other state variables have monotonic profiles. The main difference between our model and the standard one is the reversibility of the chemical reaction and the MT assumption. They are both properly taken into account in the reaction rate (2.21), and its subsequent dimensionless forms. In contrast to the standard model, the non-monotonic profiles in our case appear in all the state variables and for small values of the activation energy f . Therefore, it seems that reversibility and thermal non-equilibrium (Θ = 0) play the crucial role in generating non-monotonic profiles through their appearance in the reaction rate. It must be emphasized that the scope of this remark is restricted to ZND detonation waves, and that we do not tend to cover processes like Rayleigh flow, in which non-monotonic temperature profiles occur as well.
5.2. Case 2: Influence of f on the profile The influence of the overdrive degree f on the continuous reacting profiles is analysed for the following values of the parameters: = 1.0, f ∈ {1.08, 1.5, 3.0}, f = 2.0, (5.4a−c) and the results are presented in figure 4. As expected, different values of f lead to different terminal equilibrium states. However, there is no qualitative difference between the profiles -only non-monotonic profiles appear. This leads to the conclusion that the major influence on the qualitative structure of the profile still comes from the activation energy f , rather than from the overdrive degree f . This is a consequence of the fact that the activation energy is strictly related to the temperature sensitivity of the chemical reaction.
In the present analysis we took the value f = 1.08 with the aim to mimic the behaviour of the Chapman-Jouguet solution. It can be observed in the temperature profile illustrated in figure 4 that there occurs a temperature drop in the profile. A similar result has been obtained in (Marques Junior et al. 2015), where it was attributed to the prevalence of backward reaction (which is endothermic) for f = 1.0. However, this phenomenon appears also for f > 1.0 and yet has to be analysed carefully. . This means that the influence of a low level of activation energy opens the door to non-monotonic profiles, without any regard to the heat of reaction.
However, figure 6 reveals that, for f = 10.0, the value = 1.0 yields monotonic MT profiles, whereas higher values of yield non-monotonic profiles. Therefore, when the activation energy is high, the reaction has to be sufficiently violent in order to produce non-monotonic profiles. Moreover, by a continuity argument we may expect that there exists a threshold value of which separates monotonic MT profiles from non-monotonic ones and this aspect will be analysed in the sequel.

General remark on overdriven detonations
In view of the analysis related to the Hugoniot diagram (figure 1), and the detailed explanation of pathological and overdriven detonations, it may be observed that non-monotonic profiles of overdriven detonations exactly follow the scenario given by Khokhlov (1989) and Sharpe (1999). They pass through the point in which (4.20a) is satisfied, reach the turning point, and then converge to the equilibrium state. This behaviour is predicted in § 4 in the qualitative analysis of non-monotonic profiles. As already mentioned in the analysis of case 1, in the ST ZND model with reversible reaction, studied by Marques Junior et al. (2015), this behaviour cannot be obtained, but only the monotonic one. Therefore, the MT assumption appears to be a crucial modelling assumption for capturing the non-monotonic behaviour. In particular, for a certain range of values of the parameters, the single reversible reaction (1.1) in conjunction with the MT assumption yield the same qualitative behaviour in the case of overdriven detonation as the two-step reactions A → B → C with an endothermic phase following the exothermic one. It is also interesting to relate the turning points in the state space and the corresponding heat release. In the MT ZND model, taking into account (2.18b), the heat release is Q = z . The turning point in overdriven detonations corresponds to the maximum heat release (see Sharpe 1999  and this impliesτ p = 0 (see (4.10)-(4.14)), which corresponds to the local extremum in our study. However, there remains a qualitative difference between the overdriven detonations in the MT ZND detonation model with a single reversible reaction, and the same detonations in the two-step model -in the former case there appear more than one local extremum in the profile. This will remain an open question for future studies.
5.5. Qualitative analysis of the reaction rate The analysis in case 3 to a certain extent confirms the observation already given in case 1 -the major influence on the qualitative behaviour of the continuous reacting profile comes from the activation energy f . However, the heat of reaction becomes important, and actually decisive, when the activation energy is large and suppresses the non-monotonic behaviour.
These findings can be supported by the qualitative analysis of the reaction rateτ MT p already given in (4.17) bỹ For f small enough, the reaction rateτ MT p is large enough to induce the crossing condition and appearance of a local extremum, i.e. the appearance of non-monotonic profiles. When f is increased, the exponential factor inτ MT p is dominant and suppresses the non-monotonic behaviour. However, this may be changed if some other parameter, like , becomes sufficiently large to prevail and re-establish the non-monotonic profiles. That is what exactly happened in case 3.

Threshold in the parameter space
The main difference between MT and ST ZND detonation waves lies in the fact that in the MT case there appear two qualitatively different configurations of continuous reacting profiles -non-monotonic and monotonic ones. The non-monotonic behaviour of state variables is directly related to the MT assumption, since both crossing condition and local extremum are consequences of thermal non-equilibrium, Θ = 0. However, the parametric analysis revealed that the activation energy f has a major influence on the qualitative behaviour. Therefore, we performed a thorough numerical analysis of the profiles with the aim to determine the threshold in parameter space separating the non-monotonic from the monotonic ones. Figure 7 presents the threshold in the ( f , f )-plane for the fixed value of the reaction heat, = 1.0. It reveals that a threshold exists separating two classes of profiles. It may be observed that the threshold value of the activation energy f increases with the increase of the overdrive degree f , and this relation is approximately linear. The threshold values are the same for all the state variables ρ, u and z, but they differ for T -there occurs a temperature drop phenomenon for small values of the overdrive degree f 1.0 (up to f ≈ 1.1). This is indicated in the threshold graph by a narrow dark grey region at the bottom that corresponds to the non-monotonic profiles, but it is related to a temperature drop which occurs in the ST case as well. For other state variables ρ, u and z, this region should be regarded as light grey. Figure 8 presents the threshold in the ( , f )-plane for the fixed value of the activation energy, f = 10.0. Numerical analysis showed that the threshold in this case is not linear, and that there is an upper bound for for which the monotonic profiles appear. Also, the peculiar region of solely temperature non-monotonic profiles is found like in the previous case. It has to be noted that our analysis was performed for the value of f for which two kinds of profiles exist; e.g. for f = 2.0 there exist only non-monotonic profiles in the region of parameter space in our study. Îf FIGURE 8. Threshold behaviour in the ( , f )-plane: non-monotonic profiles (white region) and monotonic profiles (light grey region) of all state variables; non-monotonic profiles of the temperature T (dark grey region). Finally, figure 9 presents the threshold in the ( f , )-plane for the fixed value of the overdrive degree, f = 1.5. Like in the first case, the threshold is approximately linear, and there is a lower bound for the activation energy f for the existence of monotonic profiles. However, since f > 1.1 there is no region in the parameter space in which only temperature non-monotonic profiles appear.

Summary and forthcoming problems
In this work, we proposed a variant of the ZND detonation model for a chemically reactive mixture with the aim to study the non-equilibrium effects in the detonation wave. The standard ZND model is generalized in two aspects: (i) following Marques Junior et al. (2015), we introduced the reversible chemical reaction, and (ii) we considered a multi-temperature (MT) model for the gaseous mixture established within the framework of extended thermodynamics with the reaction rate derived from the kinetic theory.
One of the essential elements of the ZND detonation model is the finite length reaction zone following the shock wave, where the transformation of the mixture constituents to the detonation products occurs, due to the explosive chemical reaction and consequent energy release. The non-equilibrium effects in this zone determine the structure of the detonation wave solution and the model proposed in this paper allows us to capture many of these non-equilibrium effects, in particular those of thermal nature.
In fact, the MT model of the reactive mixture, combined with a detailed secondorder rate law deduced from the kinetic theory, provides a comprehensive description of the detonation solution in terms of the underlying parameters in the detonation process, namely the overdrive degree, activation energy and chemical reaction heat. In particular, the MT model seems to be appropriate to describe the effects due to the transfer of internal energy from reactants to products of the reaction, or vice versa, and to the opposite transfer of kinetic energy. As an outcome, the one-step reversible chemical reaction together with the MT assumption provide a simple framework that properly describes steady overdriven detonations, and yields qualitatively the same results as a two-step reaction.
The main results obtained with the new model are of qualitative and quantitative nature, alike. Qualitatively, the introduction of the MT assumption in the reaction rate implied non-zero temperature difference within the profile and the possible existence of non-monotonic detonation profiles. The latter is a consequence of the more complex equilibrium condition, which, apart from the law of mass action, requires the equilibration of the constituents' temperatures. However, the monotonicity of the detonation profile strongly depends on the underlying parameters mentioned above, which required a thorough quantitative analysis. By numerical computation, it was shown that the activation energy f has the major influence on the monotonicity: lower values of the activation energy allow the appearance of non-monotonic profiles, whereas higher values suppress them. This result motivated the computation of the threshold in a pertinent parameter space, which separates monotonic from non-monotonic profiles. It appeared that the threshold value of f increases with the increase of the overdrive degree f . Furthermore, a quantitative analysis showed that the temperatures of the constituents were driven apart by the chemical reaction, supported by the large discrepancy in initial concentrations of the constituents. Even in the quasi-MT case, in which the MT part of the reaction rate is not taken into account, there appears a temperature difference within the profiles.
In our opinion, the inclusion of the MT phenomenon in the mathematical model enhances the description of the detonation profiles of the ZND type. It reveals certain substantially new aspects of the solution that were not captured with the single-temperature (ST) model. Since the MT phenomenon is the hallmark of thermal non-equilibrium in gaseous mixtures, we expect that its influence on the processes in reactive mixtures still has to be thoroughly studied. This opens the perspective for further applications of the MT model. In particular, it may be proved to be a key tool in the following problems: (i) to investigate the stability of the detonation wave and the effects of the chemistry on stability; (ii) to extend the MT approach to detonation and linear stability in a quaternary mixture with a bimolecular reversible reaction; (iii) to describe in more detail the structure of detonation profiles, in particular the transition from the induction zone to the reaction zone; (iv) to study pathological detonations, or other non-ideal configurations of detonation waves; and (v) to study the spherical and cylindrical MT ZND detonations.
These are the planned topics of our future studies on detonation and stability. Some of them represent challenging problems, and additional tools will certainly be required. However, we believe that the MT model proposed here can contribute to a comprehensive analysis of the problems and, in particular, to a good description of the effects of chemistry on the detonation structure and propagation of instabilities. For example, when multi-step chemical reactions or other complex chemical mechanisms are involved in the detonation process, the thermal effects due to the coexisting exothermic and endothermic reactions or due to other dissipative processes have a great influence on the structure of the detonation wave and the MT model can be used to capture these effects.