Combustion with source flows

Mathl.

Comput.

Pergamon

Modelling Vol. 24, No. 8, pp. 95-104, 1996 Copyright@1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved 08957177/96 $15.00 -t 0.00

PII:SO8957177(96)00142-2

Combustion

with Source Flows

J. F. CLARKE AND C. A. LOWE Theoretical Gasdynamics and Theoretical Mechanics Cranfield University Cranfield MK43 OAL, United Kingdom Abstract-The set of conservation equations that governs the behavior of a chemically-active material is what one might call naturally inhomogeneous in view of the need to incorporate algebraic information

within it about the rates of production

of individual chemical species. There are circumarray of algebraic terms, by including others sources of mass, momentum, and energy within the

stancesfor which it is expedient to augment this natural

that account for the appearance of distributed material. A brief indication is given of how such ideas can be made consistent with physically-realistic events in a two-phase flow. They are then used to construct a simple model of events that take place during ignition of solid propellants that burn in the gas-phase. The model can be used ss a test problem for computer codes; it also forms a basis on which to build more detailed models of the complex physical processes involved in two-phase combustion.

Keywords+asdynamics,

Solid-propellants, Ignition,Chemical-kinetics.

1. INTRODUCTION Two questions are addressed in the present paper. One is to do with the construction of theoretical models of ignition of solid materials that burn, and release most of their energy, in the gas-phase. The other concerns the matter of validation of the sort of large-scale computational codes that will be necessary for the acquisition of numerical data from physically and chemically accurate models of the complex processes that are involved in real, two-phase, mixtures of energetic solid materials and the gases that evolve from them. The paper starts with a brief account of the association between a model of two-phase flow and the presence of smoothly distributed sources of mass (including mass of individual reactants), momentum, and energy in an Euler formulation of gaseous behavior. Based on these ideas, a simple model of the combustion processes emerges in Section 3.2, which allows one to study ignition events in the gas. Combustion in a spatially uniform system is analyzed in Section 3.3. A certain amount of progress can be made with the aid of classical analysis, and a number of basic features of the system are revealed in Section 3.3.1 and Section 3.3.2, including the initial conditions that must be set up to ensure that burning will be self-sustained. Some solutions of the full set of model equations are described in Section 3.4 for a system that is init,iated Combustion and gasdynamical waves by an igniter at one end of the combustion chamber. coexist and interact with one another in such a situation and solutions have to be obtained by numerical/computational means. The work described in this paper has been carried out with financial support from DRA, Fort Halstead, the authors are extremely grateful.

95

for which

96

J. F. CLARKE AND C. A. LOWE

The models described here are intrinsically interesting, but can also be used as test problems for validation of the computer codes that must be used to carry the work forward into more physically and chemically demanding situations.

2. CONSERVATION

EQUATIONS

If U and Fj (j = 1, 2, 3) denote the column vectors of conserved quantities? and their fluxes2, with vectors Sch and S, to acknowledge the presence of sources in the field, the conservation equations for material confined within a volume V surrounded by a closed surface A can be written in the following economical form, where nj is a unit vector normal to the surface A, UdV=-

Fjnj dA +

Jv

(Sch + Sa) dV,

u=(g,+($;,i, s&=(g7 %=(f).(2) The symbol Sh represents a column vector that signifies the fact that chemical species can be created or consumed by chemical reactions, the instantaneous rates R, of which will generally vary with both time and location within V. It is clear that S& represents a collection of natural events. In the present paper, this set of natural sources will be augmented by the addition of another set, S,, which is taken to signify the existence of sources of mass &, momentum 8’, and energy fi, as well as of individual chemical species A&, the latter literally in addition to the naturallyoccurring effects of homogeneous chemical reactions R,. As such, one must regard the sources represented by S, as artificial, since it is generally true that no physical mechanisms exist for the spontaneous generation of these properties. Despite this fact some rational reasons for the presence of S, in some problems can be found, as indicated in the following sections. 2.1.Equations

for the Gases in a Two-phase

Flow

Conservation equations for a complete system of gas and particles can be given in the form of equation (l), except that one should omit the artificial source terms S, since one is then dealing with a proper physical system within which mass, momentum, and energy are not created. Suppose that V consists of the sum of V,, a volume occupied solely by gaseous material, and C Vn which, with C taken over all n, accounts for the total volume within V that is occupied by a large number n of solid particles whose individual volumes are given by V,. Each particle has a surface area A,, and the whole assembly of solid particles, plus the gas that surrounds them, is contained within the closed surface A. It can be assumed that V, and hence also A, does not depend upon the time t. The version of (1) that applies to just the gaseous part of the foregoing system reads as follows: d z

J

V-E v,

UdV=-

J A+E A,

Fjnj dA +

J V-E v,

&h

dV,

(3)

where V - C V, is synonymous with V,, and nj points out of this latter volume. Both positions and sizes of the solid particles will generally change with time so that, exchanging the order of ‘Notation is 88 follows: p is density, ca is meee fraction of chemical species a (a = 1,. . , s), ui (i = 1,2,3) is flow velocity vector, E E e + UjUj/2 is energy per unit mesa of material, e = e(p, v, h) is specific internal energy, 8 function of pressure p, specific volume v = p-l, end the set of c, quantities. The Cartesian tensor summation convention is used where necessary. 2Fluxea are given only for the Euler model of fluid motion; Navier-Stokes equations require the addition of diffusive terms to the Euler values.

Source Flows

97

differentiation and integration on the left-hand side of (3), this relation can be rewritten in the form

J

v-xv,

aU+-

pjnj dA + C J,

J

EJt

A

(~j - uun,) nj d-4 + J,_, n

%hdv

v n

where unj, which appears as a consequence of the time-dependence

(4)

of the V,, quantities, is the

velocity vector of a point on the surface of the nth solid particle relative to the fixed volume V. Evidently the second term on the right-hand side of (4) describes the net effect of fluxes into the gas-phase across gas/particle interfaces.

2.1.1.Averaged relations Calculations of behavior in gas/particle

mixtures require the solution of problems that are

essentially three-dimensional and unsteady in character. As such, they pose considerable difficulties, even for computational methods, and it is arguable that such highly detailed descriptions of system behavior are not always desirable, even supposing them to be achievable. It is therefore expedient to take formally exact relations such as (4) and simplify them via some suitable averaging processes. Two moderately recent and pertinent works that deal with such matters in detail are those by Gough and Zwarts [l] and by Drew [2]. Present purposes will be served by a very brief indication of how matters turn out. It will be assumed here that major spatial changes within V only take place in the single, axial, direction 21 = 2, and that as a consequence average values taken over cross-sections in planes ~2x3 normal to z provide theoretical models of adequate accuracy. This means that the interfacial-flux terms C JA, (8” - Uu,j)n,dA, in (4) must be ‘spread’ or ‘smoothed’ across a cross-section, as indeed must the influence of the chemical source terms Sch. Evidently the interfacial fluxes under these conditions will appear in the guise of sources that are additional to &hi in other words they will make up the elements of a column vector S,,. Discussion of many important matters pertaining to the theoretical model, such as ‘porosity’, or fraction of the system occupied by gas alone, have been avoided here. They have been addressed in [1,2], and they are not significant matters in the context of the present paper. The problem to be described below deals with events during ignition of a solid propellant; at early times little solid material is consumed, and a constant value for the porosity is therefore acceptable in such a case. In the circumstances it is clear that (4) can now be written in an approximate differential-equation form as

where the tilde - indicates Note that the integral (Rankine-Hugoniot shock term vectors are bounded

an average value taken over the cross-sectional area of the system. forms of the conservation equations admit discontinuous solutions waves and contact discontinuities) since all elements in the sourcefor physical reasons. Combination of this information with (5) will be

sufficient for present needs.

3. IGNITION

AND

BURNING

OF A SOLID

PROPELLANT

The ideas described in Section 2.1 and Section 2.1.1 will be used here to establish a simple model that can describe the gasification, ignition and subsequent combustion of those solid propellants that burn primarily in the gas-phase. In the sense that fresh, unburnt, reactants are generated continuously at the interface between gaseous and solid material, all three of the foregoing processes can be thought of as taking place throughout the combustion cycle. However, the situation at early times is especially interesting since there must be a period of evolution during which a transition from a condition generated by the initial stimulus to a condition of self-sustained burning is taking place.

J. F. CLARKE AND C. A. LOWE

98

The first task is to write down some elementary facts about the chemical character of the situation. Note that the tilde -, which denotes an average value of a dependent variable in (5), is dropped from now on.

3.1. Gas-Phase Reactants Consider a simple type of solid propellant material that consists of a homogeneous mixture of a fuel species F and an oxidiser species X; both species are assumed to have the same molecular weights. When released into the gas-phase these reactants take part in the exothermic reaction

F + X --f 2P, where P is the product species. If R, (a = F, X) stands for the mass-rate of production

of species (Y in units of msss per unit volume per unit time, one can say in these

circumstances that RF = Rx = R = -@cFcX,

(6)

where J-2is the chemical frequency, or reciprocal chemical time; for example fi=P(T,p)exp

$$ (

(7)

. >

The only other gas-phase species in the system is assumed to be the reaction product P. Solid propellants are frequently ignited by injection of hot gases into the system from an igniter; in the present case igniter gases must consist solely of product material P, injected into the system at a rate 7it,. Continuing with the theme of maximum-simplicity modelling, gasification of F and X will be assumed to take place at equal rates, which means that the elements tiF and &X in S, (cf. (2)) will both have the same value, written here as (l/2&. Equations (2) and (5) now show that species conservation equations can be written down as follows;

Combining these equations with the overall mass-conservation equation pt + (pu)= = 7iz+ ti, (in view of the definitions given above, &f in S, is just the sum of riz and tis) it is now easy to show that cF = cx = c,

vx, t.

(2)

Therefore, only one differential equation is required to control behavior of mass fraction c; it will be exhibited in the next section. 3.2. Full Set of Equations for the Ignition Model A complete set of conservation and other necessary equations for the simple model of solidpropellant ignition can now be written down. Together with thermal and caloric equations of state, which will be used here in their simplest forms, namely3

pv=RT,

(10)

where R is the appropriate gas constant, and

e = eth +

(CF + Q)&

=

eth + 24,

3At high pressures, it is neceeeary to replace v in (10) by v covolume.

eth =

CUT,

(11)

b, although this will not be done here; b ia called the

Source Flows

99

where eth is specific thermal energy, Q is the specific energy of formation of the reactants and C,, is the (constant) specific heat of the gas mixture at constant volume, equations (5) translate into the set . . = m+m,,

(12)

= -pR2+

gh,

(13)

= (ti + 7h,)?.i, =(ti+i,,)(

(14)

ethfiU2+p

>

+rit(Q-L).

(15)

L is the heat that must be supplied from the gas-phase to gasify unit mass of the combustible solid. All gaseous species are assumed to have the same specific thermal energies at the same temperature, which means that reactants F and X each carry energy from the solid propellant into the gas-phase at a rate given by (I/2)ti(eth + (1/2)u2 + Q) units of energy per unit volume per unit time. The incoming reactants do work on the gas phase at a rate r&u in the same units. If one recalls the assumptions about the character of igniter gases it can be seen how the right-hand side of the energy equation is built up. For simplicity it has been assumed that the sources S, add mass to the system with a velocity that is exactly equal to the local flow velocity u. 3.3.

The Spatially-Uniform

Case

Specimens of solid propellants are often burnt under laboratory conditions in small closed vessels in order to measure their combustion characteristics. The time for passage of gas-dynamical waves across the small scales for ignition and vessel can therefore be equations (12) through

volume of these vessels is short and, particularly, shorter than the timeburning of the solid material. Spatial variations of pressure within the neglected, as can those of the other primitive variables, which means that (15) in the previous section reduce to the following:

deth

dt

= 2QRc2 - vti(L

- pv).

Initiation is to be by some external agency that is not modelled here, which means that 7iz,, = 0; ti will just be ‘switched on’ at time zero. Defining the independent variable r via dr -&=vrh, and using (16), it can be seen that p = ~eexp~,

(19)

where ps is the density at time zero; therefore

Ps(expr - 1) = m(t), where m(t) = s,” mdt * is the total mass of propellant burnt per unit volume in the time-interval 0 to t. Evidently either 7 or m(t) constitutes a useful progress variable in the present situation. Using 7 in place of t in the set of equations (16)-(18) above, and defining the dimensionless auantities

J. F. CLARKE AND C. A. LOWE

100

one finds that dc - = dr

-DC2+

dt’ ;i-; = 22)~~ - E E N(B,c), where E G & - (y - l)e,

(23)

Note that 7 if the ratio of specific heats for the gas, defined here by the expression (y- 1) = R/C,; both 7 and L: are introduced here simply for convenience. The product ~7i7 has the dimensions of reciprocal time, the latter evidently indicative of the time taken to inject unit mass of propellant into unit volume of the gas-phase. Since 52 is a reciprocal chemical time, it is clear that 2) is a local Damkohler number whose value can vary considerably with time, principally as a result of the strong temperature dependence of R. It is acceptable on physical grounds to assume that 2) depends only on the value of 8 (cf. (7) and observe that, in real situations, rh will also be connected with 0 in some way that need not be made explicit here). 3.3.1.

The (c, @-plane

The pair of equations (21) and (22) is equivalent to the autonomous differential equation &? = N/D which has singular points in the (c, @-plane wherever N and D vanish simultaneously. $his will happen where 2, = D8, c = c,, with subscript s quantities defined by relations n

8

D, =

(lf,

2’

.

2c, = 1 -E,.

(24)

.3

Only the segment 0 < c < l/2, 0 < 8 of the (8,c)-plane is of physical interest; the second of equations (24) implies that E, must be restricted to the range 0 < J& < 1 in this part of the plane. The first of equations (24) will provide a unique solution for es, and hence for cs, if s E DL > 0; it will henceforth be assumed that ‘0, behaves in this (physically acceptable) way. Linearizing (21) and (22) in the neighborhood of the sole singular point at (c,, 0,) it can be shown that eigenvalues of the Jacobian matrix of N and D are real and of opposite signs, so that (cB, 0,) is a saddle point. A sketch of the integral curves in the (c, @-plane, which incorporates the foregoing information 8s well as one or two other features of the system, is given in Figure 1. The picture is qualitatively the same under all circumstances for which the general condition D’ > 0 applies. 3.3.2.

Physical

features

of the system

Combustion systems of the type modelled here start from the condition that no gaseous reactants are present until the application of some stimulus at time t = 0. It is clear from the foregoing analysis, and from Figure 1, that the initial stimulus must ensure two things (i) that 7itis given a positive value, and (ii) that the ambient gas temperature must be raised to a value such that 8 > 81, where 01 is defined by the point of intersection of the lower separatrix from (cB,0,) with the c = O-&s. It can be seen from (21) and (22) that $lc=e = -22, so that $lc=e 5 0 for e 2 0; equality occurs when B = 811 G L/(y - 1). If 8 > 011, gas temperature will rise from the time of initiation but, if 81 c 8 < 011, thermal energy will at first diminish (until such time as the particular integral-curve crosses the locus

Source Flows

101

0 Figure 1. The physically interesting part of the (c,0)-plane. Arrows on the integral curves show the direction of t increasing. Consult Section 3.3.1 for details and definitions.

N = 0) as a consequence of the initial need to supply heat from the gas to the solid in order to sustain the system at early times. Provided that a su5cient quantity of propellant is supplied (i.e., provided m is sufficiently large), any system that can be put into a state 8 > 01, c = 0 at time t = 0 will be self-sustaining or could, indeed, be said to have ignited. It is evident that there is always a maximum value for the reactant mass fraction when the system is operating from an initial state c = 0 in a self-sustaining mode. Some numerical solutions of (21) and (22) for a particularly simple situation, in which ti is given a constant value* for t 2 0, are summarized in Figure 2, which shows c, 8 and 2Dc2 versus time t. For the chosen parameter values 2Dc2 is proportional to the chemical power generated by gas-phase combustion, and it is interesting to see the way in which the burst of power at early times correlates with the very rapid rise of thermal energy and pressure within the system (note that 8 and p are closely proportional to one another for the data used in Figure 2). It is interesting to examine the limit 2) -+ 00. For physical reasons chemical power must be finite so that, with Dc2 bounded, c must approach zero as the limit is approached. Equation (21) then shows that 22%~~-+ 1; consultation of Figure 2c shows that the numerical results have been run to the point where they have included such asymptotically-fast chemical activity. Calculations of internal-ballistic behavior often use models that, from the outset, make the implicit assumption that 2Dc2 = 1. Equation (22) then reduces to

de -&

=

(7 - 1)e + (1 - L),

with (1 - C) described as the combustion energy of the propellant. In the present very simple situation (25) has elementary analytical solutions’; variations of 0 with time then follow from such solutions once the gasification-rate h(t) is known. Clearly gasification, not gas-phase combustion, is the rate-controlling process in the model implied by the form of (25). 4rir = 60 kg/m38 and L: = 0 (in this limiting case the saddle point ls at ((l/2), 0); system behavior is qualitatively similar to the csse 0 > 811); n is calculated from (7) with EA = 5OOOR,P = p x 103s-’ and p given in Pascals. 51nclusion of co-volume effects (cf. footnote 3) means that the coefficient (y - 1) in (25) must be replaced by 6) and analytical solutions are not then so elementary. 47I)l(u-

J.

102

F.

A. LOWE

CLARKE AND C.

0.03-

0.457

c

0

-

-

0.02-

0.35-

O.Ol-

0.25-

0.0 0

I

I

I

1.0

I

2.0

I

1

3.0

I t

I

4.0

0.15

0

I

I

1.0

f

I

2.0

I

I

3.0

I

t

1

4.0

Figure 2. Variations of (a) mass fraction c, (b) dimensionless thermal energy 0, and (c) chemical power (proportional to 2’Dc?) with time t in ms for combustion in a spatially-uniform system; propellant burning is ‘switched on’ at time zero. Consult Section 3.3.2 for details. 3.4.

Ignition

Coupled

with Gasdynamics

Most often solid propellants are burnt in chambers whose size is sufficiently large to make times-of-passage of gasdynamical waves across the vessel comparable with ignition times. Spatial variations, far from being negligible, as in Section 3.3, now have a central part to play in the behavior of the whole system and it is necessary to solve the full set of conservation equations presented in Section 3.2. The only way in which this can be done is through the mechanism of numerical/computational techniques, and some results that have been acquired in this way will be discussed shortly. The particular methods deployed in the present case involve a combination of time-operatorsplitting (to cope with the influence of the sources) and estimation of fluxes J” through the use of a Godunov-based scheme to account for gasdynamical actions. Details of the algorithms can be found elsewhere [3]. The present problem assumes that events take place in a vessel, of uniform cross-section and length 0.762m in the z-direction, which initially contains a uniform distribution of propellant material and an igniter that occupies the left-hand end of the tube between 0 5 z 5 0.127 m. The igniter mass flux 7iz, (= 40 kg/m3s) is switched on at time t = 0, and is held constant thereafter. Gasification of the propellant is assumed to be switched on when local temperature reaches 444 K,

Source Flows

103

and 7it is also assumed to remain constant from that moment on. Other parameters have the values that have been used to construct Figure 2 (cf. Section 3.3.2). These rather naive proposals are not hopelessly out of touch with the complicated activity that is found in real situations and will serve for present, purely illustrative, purposes. Figure 3 shows spatial distributions of pressure p, reactant mass fraction c, temperature T and chemical power (proportional to ac2) for three times after initiation of events within the fixed vessel. Pressure and temperature, the values of which are typical for situations of this kind, are given in dimensional form. Information that the igniter has been switched on is transmitted into z > 0 by a primary shock PS (cf. Figure 3a), that reflects from the right-hand end of the chamber and travels back towards 2 = 0 as reflected shock RS. For the two early times, 0.75 and l.Oms, Figure 3b shows that c is zero for some distance behind PS, but then rises rapidly in the regions around z = 0.3 m and 0.6 m, respectively, until, as can be seen in Figure 3d, homogeneous chemical reaction begins in earnest. Compare this behavior with the situation described in Section 3.3.2 and illustrated in Figures 2a and 2c. Reactants continue to appear in the gas-phase to the left of the sudden rises in c and it is apparent, from comparison of the results for times t = 0.75 and l.Oms with those for t = l.!jms, that pressures and temperatures continue to rise in the left-hand sections of the vessel without any help from reflected shock RS. The latter certainly provokes high temperatures and pressures

0.4

0.6

08

X (4 C

1 t

0.02-

1: ,’ r’

0.0 lI I I I I

00

0.2

. . . . . . . . . . . . . .. 1 0.4

’

\ ; \ \ I I I I I OG

. . . . . . . -.f I

I 0.R

X

(b) Figure 3. Spatially-nonuniform ignition of solid propellant in a cylindrical vessel. Profiles of (a) pressure p in MPa, (b) reactant mass fraction c, (c) absolute temperature T in degrees K, and (d) the product fIc2 (gss-phase reaction rate) in units of s-l; times t = 0.75ms (full line), lms (dashed line) and 1.5ms (dotted line); axial coordinate CC in m. The system is initiated by an igniter situated in 0 < 1: < 0.127 m; details can be found in Section 3.4.

J. F. CLARKE AND C. A. LOWE

104 1

. ...* . . . . . .._.......... . ,oo!~ ; :,

T

t

2000

. .. .

+---7---l

\

I

I

0.2

0.0

:

‘.

e---x

-______---

__...

\

.

2.

I

0.4

.. *.. .

--\-__

I1 0.6

0.6 X

(c)

1

nc"

,\; '\!..\! ._..,....*_ 1: I 1 \..."“'.. cc

zoo-

-r-

\ \ \

0

I-IT 00

0 04

0.2

I

I 06

**..... 8

X

1 06

(4 Figure 3. (cont.)

behind it. If emission of light correlates with high temperatures,.the results in Figure 3 might propose that such a thing will first be observed at the opposite end of the vessel from the igniter; some informal observations of events in a test-rig are consistent with such a proposal. Behavior of pressure and temperature over the domain occupied by the igniter (cf. Figures 3a and 3c) is strongly reminiscent of what is seen in the simple source-driven chemically-inert flow described in (41. The exact solution of a source-driven Euler flow, that is given in [4], is also particularly useful as a test problem for the kinds of algorithm used to compute solutions like those discussed in the present sections.

REFERENCES P.S. Gough and F.J. Zwarts, Modeling heterogeneous two-phase reacting flow, AZAA Jnl. 17, 17-25. (1979). D.A. Drew, Mathematical modeling of two-phase flow, Ann. Rev. Fluid Me&. 15, 261-91, (1983). C.A. Lowe, E.F. Toro end J.F. Clarke, Numerical problems from propellant systems, Pmceedings of the First Asian CFD Conference, Volume 2, pp. 541-548, January 1995, Hong Kong University of Science & Technology. J.F. Clarke, Compressible flow produced by distributed sources of mass; an exsct solution, CoA Report 8710, Cranfield University, (1987).