## Abstract

A competitive reaction system is considered, under which some chemical reagent decays by means of two simultaneous chemical reactions to form two separate inert products. One reaction is exothermic, and the other is endothermic. The governing equations for the model are presented, and a weakly nonlinear theory is then generated using the method of strained coordinates. Travelling-wave solutions are possible in the model, and the temperature is found to have a classical sech-squared profile. The stability of these moderate-amplitude temperature solitons is confirmed both analytically and numerically.

## 1. Introduction

Nonlinear travelling waves of permanent form in a channel of water have been the subject of intense interest for a long time. When the waves are of moderate amplitude and are moving at close to the critical speed, it is known that their shape may be described approximately by the Korteweg–de Vries equation. This famous partial differential equation is discussed in detail by Whitham [1], for example, and in numerous research monographs such as the recent work of Vanden-Broeck [2], §5.2. There is now a vast literature on this equation and variants of it, and the recent overview article by Grimshaw [3] surveys some of this work. There are also other related weakly nonlinear model equations for describing water waves, and some of these generalizations are presented by Johnson [4].

There is similarly an extensive literature on travelling waves in a variety of reaction–diffusion systems. Nonlinear reaction fronts occur throughout population and disease-spread models in biology [5] and also in auto-catalytic chemical reactions, as discussed by Gray & Scott [6]. Thermokinetic models of combustion in either solid or gaseous fuel beds also predict the existence of travelling fronts of temperature, and have been analysed by Matkowsky and Sivashinsky [7] and Weber *et al*. [8], for example. More recently, there has been some interest in competitive type combustion reactions, in which a compound decays to form products through either of two possible different reaction pathways. One of these is exothermic (releasing heat), and the other is endothermic (requiring heat input), and such systems have been analysed fairly comprehensively by Ball *et al*. [9], Hmaidi *et al*. [10] and Sharples *et al*. [11]. Recently, a weakly nonlinear asymptotic solution was presented by Forbes [12], and it permitted a solution in a phase plane.

In many combustion systems, the chemical reactions cannot be considered in isolation, but instead are coupled to fluid-flow behaviour. The text by Zel'dovich *et al*. [13] discusses the basic fluid-combustion problem, and in particular describes how the rate of change in the concentration of a chemical species is due not only to reaction dynamics but also to the convective movement of the species in the fluid. Forbes [14] used such ideas in a model of bushfire spread, and an analysis of one-dimensional fluid-combustion behaviour was presented by Forbes & Derrick [15]. They demonstrated that small conflagrations can occur as soliton-type phenomena in the temperature, but that disturbances of larger amplitude could result in shock waves, even to the extent of the complete exhaustion of the fuel.

In this paper, a competitive reaction scheme is considered, in which the reagents are components of a flowing gas. The scheme is somewhat similar to that considered by Forbes & Derrick [15] and its governing equations are outlined in §2. Here, however, a weakly nonlinear spatio-temporal approximation is developed in §3 and makes use of the method of strained coordinates. Travelling waves in this system are discussed explicitly in §4, and the stability of these solutions is analysed in §5. Some concluding observations and suggestions for further investigation are given in §6.

## 2. The combustion model

Consider some chemical species X, which decays by one of two possible mechanisms. The first is an endothermic reaction, at some temperature-sensitive rate *k*_{1}, to give an inert product A. This is in competition with a second reaction, at rate *k*_{2} and occurring exothermically, under which X forms an alternative inert product B. This reaction scheme may be represented in the form
2.1A scheme of the form (2.1) was considered by Hmaidi *et al*. [10], and travelling waves in these competitive reactions have been studied by Ball *et al*. [9] and Sharples *et al*. [11]. Recently, a weakly nonlinear asymptotic investigation of these travelling waves was presented by Forbes [12] in an analysis that permitted a complete solution in a phase plane. All these studies involved the pure chemical system (2.1) and its associated thermodynamics alone; in this present study, this is now generalized to include the effects of the reaction in a compressible gas.

The competitive reaction (2.1) occurs in a gas of total density *ρ*. For simplicity, it is assumed that the flow is essentially one-dimensional, with local gas speed *u* in the direction of the *x*-axis. The rate equation for the chemical species X in equation (2.1) then becomes
2.2in which [*X*] denotes the molar concentration of species X. For the present, it will be assumed that the temperature-sensitive reaction rates follow Arrhenius kinetics (see [6], p. 84) and have the forms
2.3Here, *E*_{1} and *E*_{2} are activation energies for the two reactions, and the constants *Z*_{1} and *Z*_{2} give the rates of reactions at some standard temperature. The quantity *R* is the universal gas constant.

Conservation of mass within the gas requires that the continuity equation
2.4be satisfied. In addition, the Euler equation, expressing the conservation of linear momentum in the absence of viscous diffusion, takes the form
2.5Here, *p* is the pressure and it is related to the temperature *T* by the gas law
2.6Finally, there is an equation expressing the conservation of energy. This may be derived from the first law of thermodynamics, as described by Liepmann & Roshko [16], p. 190. It follows that
2.7where *k*_{1} and *k*_{2} are the temperature-sensitive rates (2.3), and *Q*_{1} and *Q*_{2} are the enthalpies of exothermic and endothermic reactions per mole, respectively. The constant *γ* represents the ratio of specific heats, and for a diatomic gas its value is approximately 1.4. The first term on the right-hand side of (2.7) represents Newtonian cooling to ambient temperature *T*_{a}. This form (2.7) is consistent with that given by Anderson [17], p. 196 and Zel'dovich *et al*. [13], p. 238, and the first two terms in parentheses on the left-hand side correspond to the total enthalpy in the system.

Non-dimensional variables are now introduced, based on the time and energy scales appropriate to the exothermic reaction in the scheme (2.1). Thus, the unit of time is taken to be 1/*Z*_{2} and lengths are scaled relative to with these two parameters taken from the Arrhenius rate formula for *k*_{2}(*T*) in equation (2.3). The temperature is scaled using the quantity *E*_{2}/*R* and speeds are referenced to the term . Far ahead of the reaction front, the gas density has some value *ρ*_{0} and this is taken as the reference for density. The pressure is scaled relative to *ρ*_{0}*E*_{2}, and concentrations are made dimensionless by reference to *ρ*_{0}*E*_{2}/*Q*_{2}. There are thus five key dimensionless parameters in this problem, and these are
2.8The constants *α* and *ϵ* are respectively the ratio of rates and activation energies for the two competing reactions and *γ*_{Q} is the ratio of their enthalpies. The dimensionless rate of Newtonian cooling is *β* and *θ*_{a} is the ambient temperature. In addition to these five parameters, the solution is also dependent upon the ratio *γ*≈1.4 of specific heats for the gas.

In these new dimensionless variables, the rate law (2.2) takes the form
2.9in which the Arrhenius relations (2.3) become simply
2.10The continuity equation (2.4) and momentum equation (2.5) retain the same forms in the new variables, with the gas law (2.6) becoming simply *p*=*ρT*. Finally, the energy equation (2.7) now becomes
2.11The appropriate conditions far ahead of the reaction front are
2.12in these non-dimensional variables, and the problem is considered over the entire real *x*-axis.

## 3. Weakly nonlinear approximation

In this section, an approximate set of equations is derived, under the assumption that the variables differ from the steady-state conditions (2.12) only by some small amount, characterized by a parameter *κ* assumed to be small. In addition, the two dimensionless rates (2.10) for the competing reactions are approximated with expressions of the form
3.1These relations (3.1) mimic the high-temperature behaviour of the Arrhenius rates (2.10), but they avoid the ‘cold boundary’ problem discussed by Matkowsky & Sivashinsky [7] and Gray *et al*. [18], in which the reactions would proceed at some finite rate for any temperature above absolute zero. Similar forms to (3.1) were used by Forbes & Derrick [15].

The approximation technique is based on the method of strained coordinates, an example of which is given in the text by Jordan & Smith [19], p. 148. New time and space coordinates are defined by means of the relations
3.2in which *κ* is the small parameter and *μ* is an unknown constant.

The dimensionless equations in §2 are now written in terms of the strained coordinates (3.2). The dependent variables, along with the constant *μ* in equation (3.2), are represented by means of expansions in powers of the small parameter *κ*, according to the expressions
3.3

These expansions (3.3) are substituted into the governing equations, and terms at each order of the parameter *κ* are collected. Without loss of generality, it is possible to set *μ*_{0}=1, and this is assumed from now on.

The modified reaction rates (3.1) yield the expressions
3.4The rate law (2.9) and the energy equation (2.11) are now expanded in powers of the small parameter *κ*, as in equations (3.3). At first order in *κ*, both these equations require that
3.5It then follows from the gas law *p*=*ρT* that
3.6in view of the relation (3.5). The mass equation (2.4) and the energy equation (2.11), combined with the expressions (3.4) then show that the second-order temperature term is expressed in terms of the first-order speed function by means of the relation
3.7in which the constant
3.8has been defined for convenience.

The momentum equation (2.5), with the mass equation (2.4), and the relations (3.6) give rise to the condition
at first order in the small parameter *κ*, and in fact, a similar relationship holds for the second-order term *U*_{2} in the system (3.3). This relation suggests that a travelling-wave solution should exist for this problem, and this is examined further in §4. After some algebra, it becomes evident that the first-order speed function *U*_{1} must satisfy a compatibility condition at the second order in *κ*. This relation has the form
3.9in which *Γ*_{1} is the constant defined in (3.8). This equation (3.9) can be integrated with respect to scaled time , using the upstream boundary conditions (2.12), to give the weakly nonlinear partial differential equation
3.10for the perturbation speed . The second-order temperature function *T*_{2} can then be calculated from equation (3.7).

## 4. Travelling waves

The governing relation (3.10) in the weakly nonlinear approximation is essentially Burgers' equation (see [1], ch. 4), and it yields solutions in the form of travelling waves. These are obtained by making use of the moving coordinate
4.1in the strained system. At zeroth order in the small parameter *κ*, the speed in these coordinates is therefore .

In this moving frame (4.1), the partial differential equation (3.10) reduces to the ordinary differential equation This expression may immediately be integrated once to give 4.2using the upstream boundary conditions (2.12). The relation (4.2) now admits a solution in the final form 4.3with no loss of generality. Essentially, this same form (4.3) was found by Forbes & Derrick [15] in a related combustion problem.

Equation (3.7) for the perturbation to the temperature takes the simpler form 4.4and when this is combined with the solution (4.3), it yields the temperature profile 4.5The temperature may therefore be represented as 4.6and has the form of a classical soliton (see [1], p. 468). Similar temperature profiles were observed by Forbes & Derrick [15].

It is interesting to re-construct the governing differential equation for the second-order temperature perturbation in equation (4.5). This is done in a straightforward manner, using equations (4.2) and (4.4). After some algebra, the result is
4.7It can be verified by direct differentiation that (4.5) is a solution of this differential equation (4.7). Nevertheless, the properties of this equation (4.7) are very different from those of the famous Korteweg–de Vries equation, even although both admit the soliton-type solution (4.5). In particular, the equation (4.7) has only one equilibrium point (*T*_{2},*T*′_{2})=(0,0) in the phase plane, and it is a degenerate saddle. This may be confirmed by using the solution (4.5) to demonstrate that
at the origin in the phase plane; thus two separatrices emerge from the equilibrium point, but their slope depends on the parameter *μ*_{1} that defines the amplitude of the soliton and is arbitrary. Figure 1 shows the (*T*_{2},*T*′_{2}) phase plane for solitons of 10 different amplitudes, at a representative set of parameter values.

## 5. Stability of solitons

It is possible to assess the stability of the travelling waves in §4 directly from the temperature equation (4.7), but it turns out to be easier to work from the equation (4.2) for the velocity perturbation function *U*_{1} and this is presented here.

In the travelling-wave coordinate (4.1), stability of the solution (4.3) is assessed by adding some perturbation to create
and substituting this into the governing equation (4.2), retaining only first-order terms in . This yields the linearized variational equation
5.1This stability equation (5.1) admits a general solution in the form
5.2in which *C* is an arbitrary constant. When the solution (4.3) is used in equation (5.2) and the integral evaluated, it may be shown that
5.3Consequently, the ratio in equation (5.3) is always less than 1 and in fact decays to zero as , so that the perturbation decays to zero. Thus, the original profile (4.3) is stable.

Stability of the profile (4.3) in the more general spatio-temporal context has also been examined, using the method of lines to solve the partial differential equation (3.10) numerically (see [20], p. 302). The spatial derivatives are approximated using second-order central differences, and the resulting ordinary differential equations are integrated forward in time, using the fourth–fifth-order Runge–Kutta Fehlberg quadrature in the *MATLAB* package. The system is given the initial condition
5.4for some small parameter *q* and appropriate perturbation , and integrated forward in time. If stable, the effect of the initial perturbation *f* will diminish in time, leaving simply the travelling wave (4.3).

The results of a sample integration are shown in figure 2. For comparison, the parameters are the same as those in figure 1. The initial soliton amplitude in the condition (5.4) was *A*_{T}=0.01, and the parameter *q* was chosen to be . The perturbation function used in this figure was , and 501 spatial points were used, over the interval . The initial perturbed velocity profile is labelled on the diagram, and solutions at the 20 later times are shown on a three-dimensional plot. It is clear that the travelling-wave solution (4.3) is re-established within a very short time, confirming that the velocity profile is stable, in accordance with the linearized stability analysis presented in equation (5.3).

The second-order temperature perturbation function may be computed from the velocity profiles in figure 2, using straightforward differentiation, following equation (3.7). This is achieved here using centred differences in the method of lines solution scheme. Results for the same case as in figure 2 are displayed in figure 3, at the same times , although the temperature at the initial time is not shown. The first two profiles shown are still influenced by the effects of the initial disturbance, but it is evident that the solution quickly re-establishes the travelling soliton form (4.6). Thus later times show a sech-squared profile of permanent form, that moves with speed for this choice of parameters, and this is confirmed by a careful inspection of the results.

## 6. Conclusions

This paper has presented an analysis of combustion waves of small amplitude, for a competitive reaction system in a flowing compressible fluid. A weakly nonlinear development, based on the method of strained coordinates, shows that the temperature disturbance may propagate as a classical sech-squared type soliton, although the governing equation for temperature is not the classical Korteweg–de Vries equation. This soliton has been shown numerically and analytically to be stable, at least to the order of approximation for which this theory is valid.

Nevertheless, as with the classical Korteweg-de Vries equation, the theory is only valid for the conditions under which it has been derived, and too much should not be expected of it otherwise. In particular, Forbes & Derrick [15] showed that, when the amplitude of the disturbance becomes too large, a more complete theory indicates that the temperature soliton ceases to be stable, and instead develops into a shock wave. Their model did not include the effects either of fluid viscosity or temperature diffusion, as is also the case here, and so it remains an interesting open question as to what extent diffusion might render the travelling temperature soliton stable, even for larger-amplitude waves.

The second-order temperature perturbation function *T*_{2} propagates as a classical sech-squared soliton, although it does not satisfy the Korteweg–de Vries equation. It is possible to derive the weakly nonlinear equation for *T*_{2} by combining equations (3.10) and (3.7); this is not shown here in the interests of space, although the reduced version of it in travelling-wave coordinates is given in equation (4.7) and its phase-plane behaviour is displayed in figure 1. It is therefore an interesting question as to how solutions to the spatio-temporal equation behave more generally. One such numerical experiment is shown in figure 4, for the same parameter values as in the previous figures. Here, equation (3.10) has again been solved using the method-of-lines technique as previously, but now with an initial condition consisting of *two* disturbances of the form (5.4), but centred at and added. For illustrative purposes, *L*=20 and the amplitude of the forward disturbance at is taken to be 70 per cent of the rearward disturbance. The temperature perturbation *T*_{2} has then been computed from the method-of-lines solution, in this case with 2401 spatial mesh points over the interval , and is presented in figure 4. This diagram shows that the larger rearward soliton moves faster than the forward one and eventually catches up to it, as time progresses. However, unlike the celebrated result for the Korteweg–de Vries equation, for which interacting pulses pass through one another and regain their identity after the interaction is complete [21], in the present reaction the two temperature solitons evidently coalesce to form a single larger pulse. A more systematic study of the interaction between temperature pulses in combustion systems of this type may therefore prove worthwhile, and is left for future research.

- Received October 8, 2012.
- Accepted November 21, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.