## Abstract

In this paper, the properties and stability of combustion waves propagating in the composite solid energetic material of the shell–core type are numerically investigated within the one-dimensional diffusive-thermal model with heat losses to the surroundings. The flame speed is calculated as a function of the parameters of the model. The boundaries of stability are determined in the space of parameters by solving the linear stability problem and direct integration of the governing non-stationary equations. The results are compared with the characteristics of the combustion waves in pure solid fuel. It is demonstrated that a stable travelling combustion wave solution can exist for the parameters of the model for which the flame front propagation is unstable in pure solid fuel and it can propagate several times faster even in the presence of significant heat losses.

## 1. Introduction

Investigation of the influence of inert elements on the burning rate of composite solid energetic materials has a long history and has been considered in several applied aspects: (i) development of solid propellants [1–5]; (ii) modification of the solid fuels [6–8]; (iii) self-propagating high temperature synthesis and the chemical furnace technologies [9–12]; (iv) development of nanostructured energetic materials and thermopower waves [13–15]. Although, the planar geometry is sometimes employed [9,12], in most cases cylindrical geometry is considered, where a cylindrical inert slab with high thermal conductivity is placed in a combustible solid fuel mixture. The basic model of such a composite has the shell–core structure i.e. inert core element is covered with the solid fuel annulus. The metal wire [1] or carbon nanotubes [13] can serve as inert central heat conducting element. Its main role is to provide heat recuperation from hot products to fresh combustible mixture. Similarly, excess enthalpy principle is also well known in combustion of gas mixtures in micro-channels [16]. This allows one to successfully increase the burning rate and temperature of the solid fuel and stabilize the flame.

In several previous papers [17,18], the diffusive-thermal models are developed and used to demonstrate these features of the thermally guided combustion waves in shell–core composites and to identify the parameter values, where flame acceleration and stabilization is possible. However, these models are adiabatic and do not include the heat exchange with the surrounding media, which is inevitable in any experimental realization of the combustion of the solid composite material.

Here, we focus on the analysis of how the heat losses affect the enhancement of the burning rate and the limits of the stable combustion wave propagation in the model of solid shell–core composite energetic material.

## 2. Model

The governing equations describing the conservation of energy and mass of fuel take the form
*ρ*_{1,2} is the density, *c*_{1,2} is the specific heat, *T*_{1,2} is the temperature, *Y* is the mass fraction of fuel, λ_{1,2} is thermal conductivity, *S*_{1,2} is the cross-section area of corresponding elements of the composite structure; *P*_{in},*P*_{out} are the perimeters of the interfacial and outer orthogonal cross section of the composite; *K* and *H* is an effective heat-exchange coefficient between two solids and fuel-ambient media, respectively; *Q* is the heat of the reaction; *T*_{0} is the ambient temperature; *E* is the activation energy of the reaction and *U*_{f}, is the flame speed.

The parameters *H* and *K* are introduced here as phenomenological ones. In the case when composite material is surrounded by a gas media, the parameter *H* can be written in the same spirit as in [19] as *H*=*Nuλ*_{g}/*R*_{1}, where *Nu* is the Nusselt number, λ_{g} is the thermal conductivity of gas and *R*_{1} is the external radius of the fuel layer. If we assume that λ_{2}≫λ_{1} and the thermal boundary and contact resistance between the solid core and shell elements of the composite are negligible, then parameter *K* can be estimated as λ_{1}/*L*_{1}, where *L*_{1} is the characteristic length of the shell in the transverse direction.

Introducing the mass fraction of fuel normalized with respect to the upstream value, *Y* _{0}, non-dimensional temperatures, *θ*_{1,2}=(*T*_{1,2}−*T*_{0})/(*T*_{a}−*T*_{0}), where *T*_{a}=*T*_{0}+*QY* _{0}/*c*_{1} denotes the adiabatic flame temperature of pure fuel; the characteristic time and length scales as
*α*_{1,2}=λ_{1,2}/*ρ*_{1,2}*c*_{1,2} is the thermal diffusivity of fuel (‘1’) and inert material (‘2’), *β*=*γE*/*RT*_{a} is the Zel’dovich number and *γ*=(*T*_{a}−*T*_{0})/*T*_{a} is the heat release parameter, the non-dimensional governing equations can be written as
*u*_{f}=*t*_{c}*U*_{f}/*l*_{c} representing the dimensionless flame velocity and the dimensionless reaction rate *ω* given by

In the system of equations (2.5)–(2.7), the non-dimensional parameters are introduced as

The boundary conditions are defined so as to correspond to the combustion wave propagating in the positive *x*-direction as

## 3. Results

The analysis of the combustion waves in the shell–core composite material is based on numerical methods. The details can be found in our previous papers, while here we briefly outline the methodology used. In order to investigate the dependence of the flame speed on the parameters of the problem, we reduce the system of the governing equations (2.5)–(2.7) to the set of ordinary differential equation by dropping the terms with a time derivative i.e. we seek the solution **u**(*x*)=(*θ*_{1},*θ*_{2},*Y*)^{T} of equations
*D*_{11}=1 and *D*_{22}=*α*, and the reaction terms are defined as **N(u)**=[*ω*−*ξ*⋅(*θ*_{1}−*θ*_{2})−*hθ*_{1}, *s*⋅*ξ*⋅(*θ*_{1}−*θ*_{2}),−*ω*]^{T} in the form of the travelling wave propagating with velocity *u*_{f} towards the negative *x*-axis direction. The boundary conditions remain intact and are defined by equations (2.10). The system of equations (3.1) subject to the boundary conditions (2.10) is solved numerically using the standard shooting-relaxation algorithms [20], where shooting is realized via the fourth-order Runge–Kutta scheme and in relaxation the two-point boundary value problem is formulated as the system of finite-difference equations on a uniform grid, which is solved with a multidimensional Newton algorithm.

The solutions obtained in this way are used to linearize the governing equations (2.5)–(2.7) near them and formulate the linear stability problem i.e. we seek the solution of the form **u**(*x*,*t*)= **u**(*x*) represents the travelling combustion wave, terms proportional to the small parameter *ϵ* are the linear perturbation terms, λ is the increment of exponential growth of perturbation. Substituting this expansion into (2.5)–(2.7), leaving terms proportional to the first order of *ϵ* only we obtain
**u**(*x*). The eigenvalue problem (3.2) is complimented with the boundary conditions

The non-stationary solutions are obtained by direct integration of the partial differential equations (2.5)–(2.7), which are solved by two methods. In the first approach, the finite-difference method of splitting with respect to the physical processes is used. Initially, we solve the set of ordinary differential equations which describe temperature and fuel concentration variations due to the reaction by using the fourth-order Runge–Kutta algorithm. As a next step, equations of heat and mass transfer are solved with the Crank–Nicholson method of the second-order approximation in space and time. As the second approach, the finite-element package [23] with adaptive time step and space mesh size is also employed to obtain the time-dependent solutions of the governing system of equations. The results obtained with both methods are compared.

The values of characteristics calculated with different numerical methods such as flame speed, critical conditions for the loss of stability, frequency of oscillations, etc. were compared and the relative discrepancy was found to be less than 0.5%. The same order of disparity in the calculated values was found when the time/space integration step was varied.

### (a) Flame speed in the non-adiabatic case

The main focus of the current analysis is in qualitative understanding of the flame dynamics in composite reactive materials rather then in modelling of some specific system, thus we discuss representative results obtained for typical values of the problem parameters. The choice of parameter values is discussed in our previous papers [17,18], where they are estimated for typical solid combustible materials and inert media. Here, we use the values within the ranges outlined in these papers, namely, we consider *α*=100, *β*=8, *γ*=0.7, *s*=0.1−100, *ξ*=0.01−100.

Several points should be made before the description of main results obtained in this current investigation. In the system containing solid fuel only in the adiabatic case the travelling combustion wave is unstable for *β*=8 [22,24–26] due to the onset of diffusive-thermal instabilities [27]. This situation is quite common for solid energetic materials [28].

In composite shell–core material, it is possible to shift the onset of pulsating instabilities to the values of the Zel’dovich number larger than 8 if the parameters of the composite structure, *ξ* and *s*, are optimized [17,18] due to the heat recuperation through the inert core. For the chosen set of parameters, the maximum flame speed is reached for *ξ*∼3 and *s*∼4. This regime of strong recuperation [18] is characterized by a substantial super-adiabatic temperature of the solid fuel in the reaction zone. Once the parameters are detuned from the optimum, for example, if *s* or *ξ* are increased, the regime of strong recuperation disappears, the temperature profiles of shell and core approach each other and the flame speed is strongly reduced. In the opposite limit, if *s* or *ξ* becomes small, the shell acts effectively as a thermostat or the role of shell becomes negligible, so that the behaviour of the system tends to non-adiabatic or to adiabatic solid fuel combustion, respectively.

The inclusion of the heat loss from the outer surface of the composite adds one more critical phenomena to the system i.e. thermal flame quenching [29]. It is also known to reduce the flame speed and to enhance the onset of diffusive-thermal instabilities [26,30] for combustion of pure solid fuels. It can be expected that the behaviour of the combustion waves in composite material in the non-adiabatic case can be substantially affected by the heat losses to the surroundings especially for the parameter values for which the heat recuperation is weak.

In order to illustrate the influence of the heat loss on the flame speed in the composite system, we prepared three cross sections in the plane of parameters *ξ* versus *s* with the lines *s*=35, 4, 0.4 and calculated the dependence *u*_{f}(*ξ*) for various *h*. The results are presented in figures 1–3. The parameters are taken so as to cover different regimes found in the adiabatic system [18] and the cases of large and small values of *ξ* and *s*.

The dependence of the flame speed on *ξ* for *s*=35, *h*=10^{−3}, 10^{−2} and 0.025 is plotted in figure 1 with the solid, dashed and dotted lines, respectively. The squares denote the stability boundaries and these data are discussed in the next section. It is seen that the flame speed is only moderately affected by the heat loss for *u*_{f} with the increase of *h*. The behaviour of *u*_{f}(*ξ*) changes qualitatively as *h* is further increased as it is shown in figure 1 with the dotted line for *h*=0.025. The function *u*_{f}(*ξ*) becomes double-valued and isola-type with two turning points for large and small values of *ξ* (the latter is not shown in the figure). There are two travelling wave solutions in this case and they exist only within a certain range of parameter values. As the heat loss is further increased, the size of the isola in the parameter plane decreases and at a certain critical value of *h* it collapses and travelling wave solutions cease to exist for any *ξ* and *s* values. It is remarkable that the travelling combustion wave can exist in the composite material for very large values of the heat loss parameter, exceeding the flame thermal quenching limit for pure solid fuel, which is found to be about *h*_{c}≈0.018 for the same choice of parameters.

Qualitatively, a similar picture is observed for the case *s*=4, which is close to the optimal value for the flame speed acceleration in the adiabatic model. Quantitatively, the isola is formed here for a smaller value of *h*=0.01 and a substantially larger flame speed is achieved. It is remarkable that the maximum of *u*_{f}≈3.1 is high even for significant heat losses *h*=0.025 and heat recuperation mechanism is effective and capable of sustaining the fast travelling wave propagation in this case. There is a notable drift of the coordinate of maximum of the flame speed towards smaller values of *ξ*: for *h*=10^{−3} maximum *u*_{f} is found for *ξ*≈2.84, while for *h*=0.025—for *ξ*≈1.66.

As it is discussed above, in the case of small values of *s* the system tends to the model of non-adiabatic solid fuel combustion, where *ξ* plays the role of the heat loss parameter. As a result the flame speed reduces significantly and the turning point occurs even in the adiabatic model [18]. In figure 3, the dependence of *u*_{f} on *ξ* is plotted for *s*=0.4 and different values of parameter *h*=10^{−3} and 10^{−2}. In contrast with the cases *s*=4 and 35 considered above, the flame speed is strongly affected by the heat loss. There is always a turning point condition as *ξ* is increased, where two solution branches representing the travelling wave solutions meet, and there are no travelling wave solutions for larger values of *ξ*. The value of the parameter *ξ* at the turning point shifts to smaller values as *h* is increased.

### (b) Stability analysis

The results of the linear stability analysis are presented in figures 1–3. It is seen that for large and moderate values of *s*, when the heat recuperation mechanisms effectively enhance the burning rate and temperature, there is a range of parameter values corresponding to the stable travelling combustion waves. These parameters are located between the full and empty squares, which show the neutral stability boundaries for the onset of pulsating instabilities. As the heat loss coefficient is increased the region of parameters, for which the combustion waves are stable, becomes smaller and at certain values of *h* the travelling wave solutions lose stability as shown in figure 1 with dotted curve. It follows from figure 1 that the maximization of the flame speed does not necessarily provide the stability, because the

More detailed analysis of the linear stability of combustion waves is summarized in figure 4, where the neutral stability boundaries are plotted in the *ξ* versus *s* plane of parameters for *α*=100, *β*=8, *γ*=0.7 and different values of the heat loss coefficient, *h*=0÷0.042. The stable parameter values in each case correspond to the points of the plane located inside the closed contours representing the neutral stability boundaries. The area encircled by these curves shrinks substantially as the heat losses are intensified and at a certain value of *h*>0.042 the contour collapses. This is the point when the travelling combustion waves become completely unstable. In the case of adiabatic combustion of pure solid fuel, the flame front becomes completely unstable against the onset pulsations at *β*≈7.2 [26]. The inclusion of heat losses shifts the critical conditions to much smaller values of *β* [26,31,32]. As it is seen in figure 4 for the composite material it is possible to sustain the stable combustion wave propagation for *β*=8 even in the presence of significant heat losses.

It is remarkable that the neutral stability boundaries have one interesting particularity, which is shown in figure 5, where an enlarged fragment of diagram presented in figure 4 is shown for *h*=10^{−3} in the *ξ* versus *s* plane. The axes are plotted here in the linear scale, while the other parameters are taken as in figure 4. The region of parameters corresponding to the linearly stable travelling wave solutions is shaded. It is seen that the stability boundary is constituted of two branches, which intersect non-tangentially forming a point where the critical curve for the onset of pulsations is not smooth.

The nature of this phenomenon is clarified in figure 6, where the real and imaginary parts of the eigenvalue of the linear stability problem are plotted as functions of *s* for *ξ*=0.7 (figure 6*a*), 0.88 (figure 6*b*) and 0.95 (figure 6*c*). The choice of *ξ* corresponds to the vertical edges of the coordinate box and the dotted line in figure 5. The other parameter values are taken as in figure 5. For *ξ*=0.7, the stability is lost as parameter *s* is decreased below the critical value 9.469 and a pair of complex conjugate eigenvalues with Reλ>0 emerges. In figure 6*b*, *ξ*=0.88, which is sufficiently close to the point of singularity of the critical curve. As in the previous case a pair of the complex conjugate eigenvalues crosses the imaginary axis at *s*=1.07 and moves to the right half of the complex plane causing the onset of pulsating instabilities. However, as *s* is further decreased this pair crosses the imaginary axis one more time at *s*=1.048 and shifts back to Reλ<0 region. Prior to this event, one more pair of complex conjugate eigenvalues crosses the imaginary axis at *s*=1.052 and moves to Reλ>0 part of the complex plane causing the instability of the travelling combustion wave. Thus, in this case there are three Hopf bifurcations. The first bifurcation (*s*=1.07) is responsible for the loss of stability of the combustion wave, while the other two bifurcations (*s*=1.048 and *s*=1.052) causes the exchange of dominating instability mode. The coordinates of these bifurcations in the (*ξ*, *s*) plane of parameters are also plotted in figure 5 with the dashed line and they form a secondary loop. The travelling wave solution is also unstable inside of this loop. For *ξ*=0.95 shown in figure 6*c*, there is a single pair of complex conjugate eigenvalues which causes the transition to instability at *s*≤1.104. However, this pair is different from the one that is plotted in figure 6*a* and is a continuation of the eigenvalues which emerge in figure 6*b* at *s*=1.052. As *s* is decreased to 1.087, the eigenvalues merge so that Im λ vanish and they become purely real for *s*<1.087. This behaviour is known for non-adiabatic premixed flames [31].

### (c) The onset of flame pulsations

At the singular point of the neutral stability curve, there is a jump in the characteristics of the instability as it is demonstrated in the previous section. Most significantly comparing figures 6*a*–*c* we can conclude that there should be a jump in the frequency of small oscillations at this point. Here, we investigate this issue in more detail. The results are presented in figure 7, where the Hopf frequency is plotted against *ξ* for *α*=100, *β*=8, *γ*=0.7, and various values of heat loss coefficient *h*=10^{−3} and 0.042. For *h*=10^{−3}, the point of singularity of the neutral stability curve corresponds to the discontinuity of the *ω*(*ξ*) dependence. The frequency of small oscillations obtained by direct integration of the non-stationary equations (2.5)–(2.7) is plotted with the ‘•’ and ‘+’ symbols. Once the neutral stability boundary is crossed in the space of parameters two different scenarios can be realized. Either stable limit cycle is formed, which is characterized by the oscillations with constant amplitude, or the unstable limit cycle emerges, for which the amplitude of oscillations grows until the flame quenching occurs. In the first case (supercritical Hopf bifurcation) the symbols ‘•’ and the solid line are used, whereas in the latter case (subcritical Hopf bifurcation) the scenario is denoted with the crosses and the dotted line. The situation qualitatively changes as the parameter *h* is increased. For *h*=0.042, the point of singularity disappears and the dependence *ω*(*ξ*) becomes smooth. Crossing the neutral stability boundary in this case results in the formation of stable limit cycle. It is seen that the prediction of the frequency of small oscillations obtained with different approaches (linear stability analysis and direct numerical integration) are in good agreement.

The characteristic time histories of the maximum temperature of the shell for the case of stable (*ξ*=3.960, *s*=5.0) and unstable (*ξ*=3.60, *s*=4) limit cycles are presented in figures 8*a*,*b*, respectively. For figure 8*a*, the amplitude of pulsations converges to constant value after some transient initial period of time and remains unchanged for many hundreds of periods of oscillations. In the case of unstable limit cycle in figure 8*b*, the amplitude of pulsation is growing exponentially until the flame quenching occurs.

In figure 9, the critical parameter values for the onset of pulsations are plotted in the *ξ* versus *s* plane for *α*=100, *β*=8, *γ*=0.7 and different values of the heat loss parameter. The critical conditions for the onset of stable and unstable limit cycle are shown with the solid and dashed lines. The full circles and crosses demonstrate the data calculated with the direct numerical integration of the equations (2.5)–(2.7).

## 4. Conclusion

The diffusive-thermal model of combustion wave propagation in the composite solid energetic material of the shell–core structure is considered. The model is one-dimensional and includes the heat loss to the surroundings. The analysis is based on numerical methods. The main conclusions of the work are as follows.

For optimal configuration, the propagation of stable travelling combustion waves is possible for the parameter values, for which the flame fronts are unstable or do not exist for pure solid fuel, even in the case of significant heat losses. Thus the flammability limits of solid fuel combustion can be significantly extended by using the effective heat-recuperation scheme via the inert heat-conducting elements. Moreover, in the case of the composite material the flame temperature is substantially superadiabatic and the flame speed can be several times higher when compared with the corresponding characteristics of the combustion of solid fuel only.

Two scenarios of the onset of pulsating instabilities are determined. Once the neutral stability boundary is crossed in the space of parameters of the problem either stable or unstable limit cycles are formed leading to the emergence of oscillations with finite constant amplitude or to the oscillations with growing amplitude causing flame quenching. The parameter values for the occurrence of these dynamical regimes are found.

A new effect of the discontinuity of the Hopf frequency along the critical curve for the onset of instabilities is found and quantified. This can lead to interesting and complex dynamical behaviour due to the nonlinear interaction of instabilities with different frequencies near the singular point of the critical curve, which we plan to investigate in future work.

## Authors' contributions

All authors participated in formulation of the mathematical model and wrote the paper. V.V.G. did most of the numerical simulations. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

V.V.G. acknowledges the financial support from RFBF grant number 16-03-00758.

## Acknowledgements

The authors are thankful to R. V. Fursenko for fruitful discussions.

- Received December 23, 2016.
- Accepted February 24, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.