- Email: [email protected]
An integral model of turbulent fire plumes
Eighteenth Symposium (International) on Combustion
AN INTEGRAL
MODEL
The Combustion Institute, 1981
OF
TURBULENT
FIRE
PLUMES
FRANCESCO TAMANINI Factory Mutual Research Corporation, Norwood, Massachusetts 02062
The paper presents a model of fire plumes in which oxygen and fuel are allowed to be present in the average at a given height, even though chemical processes are assumed to be infinitely fast. The treatment of the composition fluctuations is derived from the k-~-g model of turbulence, in its adaptation to the case in which the variation of average properties across the plume is neglected (top-hat profiles). Predictions from the simple, one-dimensional integral model for an intermediate-scale propane flame show that important quantities, such as flame height and burning rate distribution in the flame, can be calculated w i t h accuracy close to that provided by the two-dimensional k-~-g approach. While final model validation still requires comparison with experimental measurements, which are not available at present, the model offers progress toward economical predictive techniques for simple combustion problems.
1. I n t r o d u c t i o n
Advanced models of turbulence and associated numerical techniques, ~ while justified for the purpose of improving our understanding of the details of turbulent flows, are probably not the most convenient tool for engineering calculations involving practical flows of simple geometry. Integral techniques have the potential for offering an alternate approach in which the important overall features of the flow are modeled realistically, at a significant saving in complexity of the formulation and in the time required to obtain solutions. In the integral approach such reduction in model complexity is achieved by assigning a prescribed shape to the cross-stream variation of flow properties; as a result, the model contains only ordinary differential equations, with longitudinal distance as the independent variable. Steward's work z on the modeling of turbulent fire plumes is an often referenced application of the integral technique to a situation in which both buoyancy and chemical reaction control flow development. The most serious limitation of Steward's formulation, as well as of others conceptually similar to it, 3-5 is the procedure used to predict the rate of turbulent combustion: instantaneous reaction with fuel of all 2'3 or of a constant fraction 4 of the entrained oxygen, or assumption of a stationary laminar flame front inside the turbulent plume. ~ The feature common to these three models of
turbulent combustion is the neglect of unmixedness. This is a phenomenon observed in turbulent reacting flows, where, due to the eddy-like structure of the flow, lumps of unburnt fuel appear interleaved with lumps of unreacted oxidizer. The overall rate of combustion is in this case determined by the reaction occurring at the periphery of the eddies and by the rate at which the total amount of available interface surface increases because of eddy breakdown caused by turbulence. This paper presents a proposal for improving the combustion-related aspects of the models in Refs 2-5: the proposal is an adaptation to the integral approach of the combustion model used in earlier work ~ with the k-~-g technique. The model will be discussed in its application to the case of axisymmetric fire plumes. In order to simplify the formulation, we will neglect the variations across the fire plume of time-averaged values of the flow properties (tophat profile assumption).
1081
2. M e a n F l o w E q u a t i o n s
2.1 Overall Continuity The rate of change of total mass flow, rh, with height, x, is given by: dm --
dx
= rn o
(1)
1082
T U R B U L E N T COMBUSTION
where rn o is the mass of ambient air entrained per unit height by the fire plume. Typically the entrained flow is calculated using an entrainment expression such as that proposed by Ricou and Spalding, ~ which reads: m' o - K(poW) '/2
(2)
where W is the vertical momentum of the fire plume at the height in question, and Po is ambient density. Equation (2) is conceptually the same as Morton's entrainment expression, 7 but for the added advantage that K is independent of the assumed velocity profile across the flow. Since we will refer to Eq. (2) later, we note that the entrainment coefficient K has been found to vary from .28 for forced jets ~ to .56 for thermal plumes. ~ Rather than relying on the approach represented by Eq. (2) we have chosen to calculate the entrained flow from: nio = ck~r(~ + ~,)
(3)
where I~ and p.~ are the laminar and turbulent components of viscosity and ck is a constant. The analytical solution of the laminar plume problem by Yih 9 indicates ck = 12. On the other hand, for the turbulent case the data by George et al, 8 when combined with the experimental estimate for p. made by Beuther et al, ~~ imply a value for ct close to 7. This latter value is the one used in our calculations. It should be noted that the fact that ck assumes a lower value in a turbulent flow than in the laminar counterpart has already been observed for the case of forced jets 'L'2 (ck = 5.5 and 8, respectively). As will be shown, the model allows one to estimate Ix~ and, through rt,, the intensity of the turbulent transport. Therefore, we feel that Eq. (3) can account in some fashion for the effect on entrainment due to incomplete development of the turbulence in the low part of the flame, while Eq. (2) obviously cannot. Another reason for our adoption of Eq. (3) in this work is the fact that the empirical coefficient, ck, seems to be less dependent than K on buoyancy effects: c k increases from 5.5 to 7, while K doubles in going from turbulent forced jets to thermal plumes. If Am o is the total entrained air flow and m, the initial mass flow, then the integrated form of Eq. (1) becomes: rh = m , + a m o .
(4)
2.2 Species Conservation The conservation of chemical species, i, can be written as:
d -(u dx
= rh', + u
(5)
where Y, and Y,.o are the average mass fractions in the plume and in the ambient, and m~ is the rate of production of species i per unit plume height due to chemical reaction. Species which do not participate in the chemical reaction, as in the case of nitrogen, have zero source (m', = 0); in the case of oxygen and fuel, however, m ' and th; are not zero for as long as unburnt fuel is present. Steward ~ assumed fixed stoichiometry and wrote (for Yr > 0): mo~ "' = srh; = -Yo, ohi'o
(6)
where u is the mass fraction of oxygen in the ambient and s is the oxygen/fuel stoichiometric ratio. Since initially the jet flow contains no oxygen (Yox.o = 0), Steward's assumption (Eq. (6)) implies that Yox - 0, for as long as Yr > 0 in the plume (see Eq. (5)). We will show in the following sections how we propose to allow for presence of oxygen and fuel at the same height in the plume, by accounting for the fluctuations in composition due to the turbulence. Let us now write Eq. (5) in a simpler form by introducing a Shvab-Zeldovich transformation: if oxygen and fuel combine in constant ratios, then the variable r = Yox - sYr
(7)
is a conserved quantity in that it is not created or depleted due to chemical reaction. The conservation equation for F is: d - - (rrh) = ron~;
(8)
dx
where F o = Yox.o is the ambient value of the composition parameter, F. By taking advantage of the definition of the total entrained flow used in Eq. (4), it is possible to integrate Eq. (8) for the case of uniform ambient conditions and write: F ni =
r,ni,
+ FoAm o
(9)
where F, = -sYs. ~ is the value of the composition parameter in the jet flow. For given initial and ambient conditions, Eqs. (9) and (4) can be combined to yield an analytical expression for F as a function of rh:
r = ro + ( r , - r o ) . , , / , a .
INTEGRAL MODEL OF TURBULENT FIRE PLUMES A similar expression applies to the case of the mass fraction of nitrogen (substitute Y,, Y,.~ and Y,.o for F, r~ and F o in Eq. (10)). Knowledge of the value of F is not sufficient to calculate Yo~and Ys; it is the task of the combustion model to determine the Yo~ and Ys associated with a given F. For the moment, we simply note that, once Yo~, Ys, and Y are known, the concentration of products can be obtained from: Y , . = 1 -Yo~ - Y e - Y..
(11)
since no other species need to be considered. 2.3 Energy Conservation The conservation of thermal energy requires that: d dx
(c.Tm) = q'~h+ c.Torho - q',
t12)
where T is the average temperature across the plume and cp is the specific heat (assumed constant) of the gases. The two quantities ti'~hand 0" are, respectively, the energy supplied by chemical reaction and that lost by radiative emission per unit plume height. The same procedure adopted in connection with the introduction of F, can be extended to the case of the enthalpy by defining a total enthalpy O as; O = c~T + QfYi
(13)
where Qf is the heat released by reaction of unit mass of fuel. The definition of Qf implies that: O'h = -Qsrh}
(14)
It is easy to show that the transport equation for O reads: d --
dx
(Ore) = Oom o -
(fr"
(15)
Equation (15) can be integrated to yield the following algebraic expression: Orb = O,rh~ + OoArho - &Jr
(16)
where Aqr is the total power radiated by the flame up to the height in question. Note that Eq. (16) would be identical to Eq. (9) if it wasn't for the presence of the radiation term, 5tj,: following the presentation of the combustion model, we will introduce a simple treatment of radiation which overcomes the complications associated with the difference between Eqs. (9) and (16).
1083
2.4 Momentum Equation The vertical momentum, W, of a plume in a coflowing stream of velocity U o, increases with height because of action of buoyancy: dW dx
= a,(po - p)~rb2 + ~hoUo
(17)
where a~ is the acceleration of gravity; p is the average density; and b is the radius of the fire plume. Note that the formulas introduced up to (but not including) Eq. (17) did not require a specification of the shape of the radial profiles of species, temperature, velocity, etc. inside the fire plume. A choice, however, is implicit in Eq. (17). While the model could accommodate other shapes, at the cost of an increase in complexity, we will restrict the discussion of the combustion model to the case of uniform (top-hat) profiles. In this case, the relationship for the plume radius is very simple: b = rh/X/'~pW
(lS)
3. Combustion M o d e l 3.1 The y-Plane So far we have been dealing with time-averaged quantities; having limited the scope of the analysis to such variables is the main cause for the lack of realism in Steward's model, z The behavior of turbulent flows cannot be understood unless appropriate consideration is given to the intrinsically unsteady, time-dependent nature of the flow. In other words, any turbulence model should account, in some more or less elaborate fashion, for the presence of fluctuations. Let us visualize a situation where, at a given height in the plume, the composition parameter fluctuates within a given range. The instantaneous values, ~, will be functions of time and position in the cross-section of the plume, and such that their time-space average is F. In order to evaluate the average values of species concentrations and density one needs two ingredients: 1) establishment of the relation between ~ and the corresponding instantaneous oxygen, fuel fractions (Yox, Ys), density (~) and temperature (/'); 2) determination of the probability of occurrence of all possible ~'s. If one assumes that chemical processes are infinitely fast, then fuel and oxygen cannot coexist (at the same time and at the same point); the relationship between ~ (or Yl) and z/ is easily obtained by remembering that by definition it is: r = Yox - s~s
(19)
1084
TURBULENT COMBUSTION
and that non-zero values of gox and #s correspond, respectively, to ~ greater and less than zero. With regard to the relationship between 5 and the instantaneous inert mass fraction, ~,, we note that the two normalized average variables (F F o ) / ( r , - ro) and (Y, - Y o)/(Y,.~ - Y,.o) obey the same differential equation and have identical boundary conditions. From this it follows that, for the two normalized variables, instantaneous values must also be identical to one another; the sought algebraic relation between ~ and #, is readily obtained. The concepts introduced here have already been used extensively in the literature: this is the reason for presenting them in cursory fashion. More details are available elsewhere. Determining the relationship between instantaneous temperature and ~ poses conceptual and practical problems, because of the presence of a term in Eq. (15), the radiated power q',, which has no counterpart in Eq. (8). If it is possible to assume that the radiated power is a constant fraction, • of the power produced by chemical reaction: q', = Xq'ch = - • Qcrh),
(20)
then one can define the total enthalpy as: O* = c T + (1 - •
(21)
and happily note that the term q'r no longer appears in the equation for O'. By invoking the same argument used in connection with ~,, it is then easy to obtain the algebraic relation linking t to Unfortunately, in the general case, where Eq. (20) may not be an accurate expression for t~, a unique relationship between t and ~, may not even exist. While a rigorous solution of this problem is beyond the capabilities of the approach described above, approximate solutions can be devised and have been used. la We are discussing this issue because we believe that it has implications in terms of the extent to which the modeling procedure can reliably accommodate more elaborate radiation treatments. Having made the point we now set it aside and assume, for the sake of the following application of the model, that Eq. (20) does indeed offer an accurate representation of the losses due to radiation. We will therefore calculate O ~ from Eq. (16) in which we have set Atlr = 0 and O = O ~ Once Ys is available from the statistics of the fluctuations, T can be obtained from Eq. (21). 3.2 Probability Distribution and Average
Quantities As stated previously, the evaluation of the mean values of concentrations and density requires
knowledge of the probability density function (PDF) of ~. The probability function H(~t) is defined in such a way that II (~) d~ is the probability of existence of concentrations in the range "~ to "~ + d~. We will follow an approach in which we assume that the PDF has a fixed shape; the width and the center of the distribution are determined from the_values of F and its mean square fluctuation, g = ~/2. The quantity g in turn is obtained from a transport equation as discussed in the following section. We choose the same form for the PDF already adopted in our earlier work, L'~ namely the shape given by two parabolae, which is reasonably close to the results from experiments and for which we have already worked out the necessary algebra. 14 The average value, q~, of a variable ~b(composition or density) can be evaluated through the integral:
i .,em~ "r
(22)
where 6(~) represents the analytical expression which relates instantaneous values of ili to those of ~,. In the general case, 6 (~) is a linear function, which is different for z/ greater or less than zero. Using Eq. (22) it is possible to calculate all average compositions, in particular Ys" Once Y~ is known, the burning rate per unit height is obtained from: d
m," = --dx (Y s,h)
(23)
3.3 Modeling of g-Equation Consider a mass flow of magnitude th and mean composition F (the composition parameter can be thought of as the mass fraction of an inert, since chemical reaction does not directly affect it). Let us further assume that the composition fluctuates i__~_ntime about its mean value, F, and define g = "t 2 as the magnitude of the mean square of the deviations from the mean. If one visualizes the flow as being an isolated stream, the only mechanism that can alter the value of g is the action of molecular diffusion processes, which will tend to reduce the intensity of the fluctuations until they approach zero. If, on the other hand, fresh fluid of uniform composition Fo is periodically added to the main stream, the level of fluctuations is going to increase for as long as F # F o. These two processes, molecular diffusion and entrainment of ambient fluid, constitute respectively the basic mechanisms of dissipation and production of g. We will examine now in some detail the practical implementation of these concepts by considering first the production of fluctuations associated with entrainment. An increment in mass flow from rn to rn + 8hi, due to entrainment of fluid of composition F o, brings
INTEGRAL MODEL OF TURBULENT FIRE PLUMES about an increase, 8g, of the level of fluctuations which can be calculated from:
1085
standard k-e-g model ~ as proposed by Spalding. ~" After some manipulation, the complete differential equation for g becomes:
rn (g + ~g) = [ g + (~r)~l -
-
d
E
. , ~-c~2 p -k- glrb 2 --(gin)= (F o - F ) 2too.
m + ~rh
(27)
dx
~m +(Fo-F-aF) ~
-
m+~m
(24)
The last term in the right-hand side of the equation represents the contribution to the overall level of fluctuations due to the fact that unit mass of flow contains the fraction ~rh/(rh + 8m), whose composition is Fo, while the mixture (after entrainment) has average composition F + ~F. The shift by 8F in the composition of the mixture, associated with entrainment of ~rh, is responsible for the second order increase (~F)2in the fluctuations of the original stream. The presence of mass ratios as weighing factors in Eq. (24) implies that g should be interpreted as a density-weighted quantity. Strictly speaking, all mean quantities introduced in the preceding sections should be thought of as Favre averages and not as Reynolds averages. However, a rigorous implementation of this fact would involve modeling of additional terms in the k- and ~- equations and, for this reason, we have not pursued it here (see 9 15 Bflger ). Equation (24) can be cast in differential form by dropping higher order terms; thus, the result for the case where all dissipation is neglected is:
dx
rho., = c~ ~r ~ ,
(28)
rather than the total entrained flow, as would be indicated by a more rigorous derivation of Eq. (27). ~ The distinction has some relevance only very low in the flame and we have made the change for reasons of convenience, to avoid numerical problems in the very early steps of the turbulent calculation. The evaluation of k, ~ and I~~is discussed next. 4. T u r b u l e n c e
Quantities
4.1 Turbulence Kinetic Energy and Its Dissipation The transport equation for the turbulence kinetic energy is given by:
d(krh) - -
d --
where k and e are the turbulence kinetic energy and its dissipation; while c,2 is a constant. Note that we have introduced in the source term the flow entrainment associated with the turbulent fluctuations:
= Pk + G~ - ~:~
(29)
dx (g,n) = ( r o - r ) ~ , n o .
(25)
In the absence of molecular diffusion, Eq. (25) is an exact mathematical representation of the mechanism governing the growth of g. It is worth noting that, in the case of a jet of uniform composition, F~, issuing in an environment of uniform composition, Fo, Eqs. (8) and (25) can be solved analytically and the solution is: g = (ro - r)(r - r,).
(26)
The level of fluctuations given by this solution is the maximum attainable and corresponds in the PDF plane to two delta functions at ~/ = Fo and zt = F, and no intermediate compositions. As has been pointed out, the rate of production of g can be calculated from an expression (Eq. (25)) which is exact. For the case of the dissipation of g, however, the situation is not as fortunate and we must resort to modeling. While recognizing that this point may be the subject of future improvements, we will for the moment simply extend to the plug flow considered here the formulation used in the
where Pk and G~ are total production of k due to shear and buoyancy effects respectively, and ~ is total dissipation. By analogy with the exact expression obtained for the term representing production of g (cf. Eq. (25)), we write: 1
Pk = - - (U - U o ) 2 mo., .,
2
(30)
where U is the average vertical velocity of the fire plume and U o is the vertical velocity of the free stream (in the case considered here it is Uo = 0). For the production term representing generation of k due to buoyancy we write: 17
Gk = agco X / ~ k ~rb2
(31)
*M. A. Delichatsios has brought to our attention the fact that Eq. (27) can be obtained by direct manipulation of the two-dimensional boundarylayer form of the model equation for g. His contribution is gratefully acknowledged.
T U R B U L E N T COMBUSTION
1086
where c o is a constant (=.15) and ~ is the root mean square value of the density fluctuations. Finally, the total dissipation is given by: E~ = p~Trb ~
(32)
For the dissipation ~, we follow standard practice L~8-2oand write:
- dx
-
{c,~P k + c,~G k - c,2Ek}
(33)
k
where c~, c~3 and c,2 are constants. We will later see that values for these constants, close to those considered to be "standard" L~8-2onamely c,~ = 1.44, c,2 = 1.82 and c,3 = .95, lead to reasonable model predictions.
TABLE I Properties of propane flame used as test case Jet velocity, u~ (m/s) Nozzle radius, b~ (m) Fuel jet temperature, T t (K) Ambient temperature, T O (K) Specific heat, cp (J/kg/K) Oxygen/fuel ratio, s (-) Heat of combustion, Qr (J/kg) Radiated fraction, • (-) Ambient oxygen mass fraction, Yo~.o (-) Jet fuel mass fraction, Yr (-)
3.251 .00635 300. 300. 1200. 3.636 4.64 ::< 1 0 7 .26 .233 1.
5. Model Predictions 4.2 T u r b u l e n t Viscosity The entrainment formula given in Eq. (3) contains the turbulent viscosity, p,. We obtain this quantity from the known expression, used by the k-~-g model, ts which reads: k2
p~ = c~ pV/-~op - -
(34)
s
where c = .09 is a constant. Note that the average density ~pop)~/2 is used instead of p. At this point the model is complete and its application to a specific situation requires only the choice of appropriate initial conditions. 4.3 Initial Conditions As we have indicated in previous applications of the k-e-g technique, L~3 the problem of how to start a turbulent flow calculation has yet to find a satisfactory solution. The reason is that the models used apply to fully developed turbulence and cannot deal with the initial growth in the transition region. However, the arbitrary nature of the procedures usually adopted to introduce turbulence in the flow should not justify incomplete documentation. In the calculations presented here, we have assumed initial values for k, e and g right at the exit of the jet nozzle as follows: k = 10-au~ g=0
(35) (36)
and for ~ the value that corresponded to equilibrium (Pk = Ek, cf. Eqs. (30) and (32)). The transport equations for k, e and g (Eqs. (29), (33) and (27)) then dictated the subsequent development of these quantities.
5.1 Constants Selection The model used in this work contains a number of adjustable constants, whose values can be chosen to improve agreement between theory and experiment. Of special concern to fire research is the variation with height in the plume of the volumetric burning rate. This quantity is important because convective and radiative heat release and, therefore, flow development are all linked to the location and the intensity of the chemical reaction. Since no direct measurements are available on the burning rate distribution in buoyancy-controlled flames, we have used as a point of reference the numerical predictions for the largest of three intermediate-scale propane flames, considered in our earlier work with the k-e-g technique. ~ The properties of the propane jet flame are listed in Table I and the model constants finally adopted are given in Table II. While the values of constants c,~(=1.44), c,2 (=1.82) and c.(=.09) are the same as those used in two- or three-dimensional applications of the k-e-g approach, changes in the other constants were made necessary by the particular nature of the integral approach as implemented here. So, for example, in
TABLE II Model constants Entrainment formula (Eq. (3)): g-equation (Eq. (27)): k-equation (Eq. (31)): e-equation (Eq. (33)):
Turbulent viscosity formula (Eq. (34)):
c k = 7.
ca2 = cp = c,1 = c,2 = cE~ =
.8 .15 1.44 1.82 .95
c~ = .09
INTEGRAL MODEL OF T U R B U L E N T FIRE PLUMES order to obtain a flame height in agreement with that given by our earlier predictions,~ (which in turn was found to agree with experiment), we had to reduce the value adopted for %2 down to .8. The values for c 2 quoted in the literature are in the 19 g 12 20 range 1.25 to 2.0 to 2.8 . For our earher flame predictions ~ we set c~z = 1.82. To counteract the effect of the resulting high level of density fluctuations on G,, the buoyancy production of turbulence kinetic energy, we have set c = .15. This choice should be seen against the fact that the true value of the correlation coefficient c , based on thermal plume dataY is close to .6-.7. The choice of c,3 = .95 was dictated by the desire to obtain above the flame tip the same entrainment coefficient (K = .56) indicated by thermal plume data. s
12
,
1087 --
L o l ~
9 t,,,.. 1/*? y\
1500
13oo
,,oo
Height, x(m) Fie. 2. Predictions of the integral model of average fuel and oxygen mass fractions, and temperature at various heights in the p r o p a n e flame.
5.2 Numerical Solution The solution of the model equations is quite straightforward. There are five ordinary differential equations (Eqs. (1), (17), (27), (29) and (33)) for rh, W, g, k and e, and four algebraic expressions (Eqs. (3), (10), (16) and (34)) for rh~,, F, O* and Additional relationships are used for the statistical evaluation of mean flow properties. The flame calculation described in the next section required less than 150 integration steps in the x-coordinate and less than 4 sec of CPU time on IBM 370/158.
Is~.
5.3 Flame Calculation A comparison between the predictions of variation of burning rate with height given by the full k-e-g
.ool2 -r c
,' ,;
9 ',,
Propane flame UI= 5.251 rn/s
0008
g:.ooo,~ m
,4
.8 Height, x(m)
1.2
1,6
Flc. 1. Calculated variation with height of the burning rate per unit plume height. The result of a calculation with the integral model is compared with the prediction from the standard k-(-g approach.
model ' and by its adaptation to the integral approach are shown in Fig. 1. As can be seen, the trends followed by the more complex formulation are faithfully reproduced by the simpler model developed in the present work. Even better agreement could have been obtained by additional model tuning (read: improved selection of constants). This further refinement didn't appear to be worth the trouble, since final validation of the model can only come from experimental data. Figure 2 shows the variation of oxygen and fuel mass fraction, temperature and unburnt fraction of jet fuel flow as a function of height. As can be seen, very near the base of the flame oxygen starts to appear in the fire plume, long before fuel is consumed. Furthermore, the average temperature of the fire plume reaches a maximum value of about 1450 K, which is well below the adiabatic flame temperature for this system (2023 K). The ratio between the flux of unburnt fuel (mYe) and its value at the nozzle (rh~) varies gradually from 1, at x = 0, to 0 at the flame tip. As can be seen in Fig. 2, 99% of the fuel is consumed at the height of 1.16 m. We are currently developing a novel experimental technique which will enable us to make a direct measurement of this quantity. Calculated values for vertical velocity, U, plume radius, b, and entrainment coefficient, K, are shown in Fig. 3 as a function of height. The vertical velocity first drops quite rapidly to about half the nozzle value, then decreases more gradually in the remaining part of the flame. The plume width increases steadily and so does the entrainment coefficient, which tends to approach asymptotically the thermal plume value (K = .56). The curve for K in Fig. 3 and that for the burning rate in Fig. 1 show an initial rapid rise, which is then followed by an inflection point. This occurs in the low part of the flame, where the turbulence undergoes rapid growth
1088
TURBULENT COMBUSTION
bU .2
K
Nomenclature
4
.155~N~
ag
Ent~
"4
b Cg2 r
acceleration of gravity, (m/s 2) fire plume radius, (m) constant in g-equation (Eq. (27)), (-) constant in entrainment formula (Eq.
(3)), (-) Cp C~l, C~2, C~3 C~ r
.4
Hegixh(m t, ) .8
1.2
i
g Gk
FIG. 3. Calculated variation with height of mean
vertical velocity, plume radius and entrainment coefficient as given by the integral model. in a region of numerical transition. We have verified to our satisfaction that this anomaly is not caused by numerical problems with the solution, but is dictated by the model equations.
k K m
m; Pk
q'~ 6. C o n c l u s i o n s
We have presented an integral model of fire plumes which uses the approach to the modeling of turbulent processes introduced by the k-e-g technique. Even though on the conceptual front our adaptation of the k-~-g technique offers few new ideas, the importance of this work lies in the fact that the integral model is shown to yield results, similar to those obtainable from the full 2-dimensional k-~-g approach, but at a fraction of the computer time (4 sec instead of 150 sec of CPU time on IBM 370/158 for. a flame calculation). This reduction in model complexity, with little sacrifice in terms of accuracy, is necessary for applications which require several flame calculations such as, for example, in problems of fire spread. Validation of the model will require comparison with experiment; as already mentioned, we are working on a facility in which we will be able to measure the burning rate distribution in the fire plume, as well as other important parameters. We can already anticipate that several aspects of the model will need improvements; the treatment of entrainment is one of them, particularly if one intends to use this formulation to calculate pool fires. The advantage of this approach lies in its simplicity, as demonstrated by the relative shortness of the computer program which implements it. The same computer program can easily be modified by external users to incorporate model expressions different from those described here.
q', Q~
specific heat, (J / kg/K) constants in e-equation (Eq. (33)), (-) constant in turbulent viscosity formula (Eq. (34)), (-) constant in formula for buoyancy-generated turbulence (Eq. (31)), (-) mean square fluctuations of composition parameter, ( - ) production of k due to buoyancy (Eq. (31)), (kg m / s 3) turbulence kinetic energy, (m2/s ~) entrainment coefficient (Eq. (2)), (-) total mass flow, (kg/s) rate of production of species i per unit height, (kg/m/s) ambient fluid entrained per unit height, (kg/m/s) production of k due to shear (Eq. (30)),
(kg m/s 3) rate of release of chemical energy per unit height, (W/m) radiated power per unit height, (W / m) heat of combustion of unit mass of fuel,
(J/kg) $
t,T U W X
y,,Y, "y,F hrho e
oxygen ] fuel stoichiometric ratio, (-) temperature, (K) mean vertical velocity, (m/s) vertical momentum, (kg m / s ~) longitudinal distance, (m) mass fraction of species i, (-) composition parameter (Eqs. (7) and (19)), (-) entrained mass flow, (kg/s) dissipation of turbulence kinetic energy, (m2/s 3) total dissipation of k (Eq. (32)), (kg
m/s ~) 0 09
total enthalpy (Eq. (13)), (J/kg) corrected total enthalpy (Eq. (21)),
p. P.T II p •
laminar viscosity, (kg/m/s) turbulent viscosity, (kg/m/s) probability density function of ~, (-) density, (kg/m 3) radiated fraction, (-)
(]/kg)
Superscripts -
~ '
indicates average value indicates instantaneous quantity indicates fluctuating component
INTEGRAL M O D E L OF TURBULENT FIRE P L U M E S
Subscripts f n ox pr 0 1
fuel diluent (nitrogen) oxygen products ambient conditions fuel jet conditions
Acknowledgments The author is grateful to John de Ris for useful discussions. The work was carried out under joint sponsorship from Factory Mutual Research Corporation and the National Bureau of Standards under Grant NBS G7-9011,
REFERENCES 1. TAMANINI,F., Combustion and Flame, 30, 85-101 (1977). 2. Steward, F., Combustion Science and Technology, 2, 203-212 (1970). 3. NIELSEN,H. J. ANDTAO, L. N., Tenth Symposium (International) on Combustion, pp. 965-972, The Combustion Institute, 1965. 4. BRZUSTOWSKI,T. A., in "Turbulent Combustion," Progress in Astronautics and Aeronautics, Vol. 58, L. A. Kennedy editor, pp. 407-430, American Institute of Aeronautics and Astronautics, 1978. 5. WILCOX, D. C., AIAA Journal, 13, 3, 381-386 (1975). 6. RIcou, F. P. AND SPALDING, D. B., J. Fluid Mechanics, 11, 21-32 (1961). 7. MORTON,B. R., Tenth Symposium (International) on Combustion, pp. 973-982, The Combustion Institute, 1965. 8. GEORGE, W. K., ALPERT, R. L. AND TAMANINI, F., Int. J. Heat Mass Transfer, 20, 1145-1154 (1977).
1089
9. YIH, C. S., Proceedings of the 1st U.S. National Congress of Applied Mechanics, E. Sternberg editor, pp. 941-947, American Society of Mechanical Engineers, New York, 1951. 10. BEUTHER, P. D., CAPP, S. P. AND GEORGE, W. K., "Momentum and Temperature Balance Measurements in an Axisymmetric Turbulent Plume," ASME paper 79-HT-42 presented at the Joint A S M E / A I C h E 18th National Heat Transfer Conference, San Diego, California, August 6-8, 1979. 11. SCHLmHTJNG,H., "Boundary Layer Theory," p. 184 and p. 609, Mcgraw-Hill, 1960. 12. SPALDING, D. B., "'Combustion and Mass Transfer," p. 207, Pergamon Press, 1979. 13. TAMAN1NI,F., Seventeenth Symposium (International) on Combustion, pp. 1075-1085, The Combustion Institute, 1979. 14. TAMANINI,F., "Numerical Model for the Prediction of Buoyancy-Controlled Turbulent Diffusion Flames (A Description of the F U E L J E T 2 Computer Program)," Factory Mutual Research Technical Report, Serial No. 22360-1, RC-B-47, October 1975. 15. BILGE8, R. W., Prog. Energy Combust. Sci., 1, 87-109 (1976). 16. SPALDING,D. B., Chemical Engineering Science, 26, 95-107 (1971). 17. TAMANINI,F., "On the Numerical Prediction of Turbulent Diffusion Flames," paper presented at the Joint Meeting of the Central and Western States Section of the Combustion Institute, Southwest Research Institute, San Antonio, Texas, April 21-22, 1975. 18. LAUNDER, B. E. AND SPALOINr D. B., "Mathematical Models of Turbulence," Academic Press, 1972. 19. GIBSON, M. M. ANt) LAUNt)Ea, B. E., ASME ]. of Heat Transfer, 98, 81-87, 1976. 20. TAMANINI, F., ASME J. of Heat Transfer, 100, 659-664, 1978.
COMMENTS S. Galant, Societe Bertin, France. In largely reducing computational times, your work allows for easy sensitivity analysis of the results of the model parameters. Have you carried out such an analysis for the k-~-g model parameters? If so, how do they affect numerical results? Author's Reply. We have not carried out a systematic analysis of the effect of changes in the values of the various constants. The discussion given in the paper on the criteria followed in selecting such
values, indicates that we have fine-tuned the model only to a limited extent. The flame height depends rather strongly on the value of ce~ and we have taken advantage of this fact to fit the integral model to the results from the 2D calculations. Any additional model optimization will have to be done by comparison with experimental data. I should point out that, in any case, this process of optimization cannot use all the degrees of freedom given by the seven (7) model constants. In fact, it would be hard to justify the introduction of drastic changes, which
1090
TURBULENT COMBUSTION
w o u l d m a k e the set o f c o n s t a n t s m u c h d i f f e r e n t from the o n e generally a c c e p t e d a m o n g k-c-g modelers.
fire p l u m e are very weakly d e p e n d e n t on the s h a p e of the P D F . T h i s c o n c l u s i o n h o l d s also for the integral m o d e l p r e s e n t e d in the paper.
H. Milter, Harvard University, USA. D i d y o u assume a s h a p e for y o u r P D F , or solve a P D F e q u a t i o n ? If the former, t h e n of c o u r s e y o u have in fact more t h a n the 7 p a r a m e t e r s y o u appear to have; y o u r results will p r e s u m a b l y d e p e n d on the a s s u m e d P D F . C a n y o u say, at this time, w h a t is the s e n s i t i v i t y of the results to the a s s u m e d s h a p e ?
Author's Reply. As i n d i c a t e d in the paper, we a s s u m e that the P D F in the 7y-plane is u n i q u e l y d e f i n e d by each pair of v a l u e s for F a n d g. O u r p r e v i o u s work w i t h t h e 2I) v e r s i o n of t h e k-e-g m o d e l (Ref. 1) h a s s h o w n that the overall features of the
1. Creighton, Lawrence Livermore Laboratories, USA. F r o m t h e oral p r e s e n t a t i o n , it w a s n o t clear to m e w h i c h c o m p u t e r code or n u m e r i c a l m e t h o d was u s e d to calculate the K-c-g model.
Author's Reply. T h e details of o u r 2 D calculations u s i n g the k-~-g m o d e l are reported in Ref. 1. T h e c o m p u t e r code a d o p t s the s o l u t i o n a l g o r i t h m for parabolic flaws d e v e l o p e d b y S p a l d i n g a n d coworkers at I m p e r i a l College a n d k n o w n u n d e r the code n a m e of G E N M I X .
AN INTEGRAL
MODEL
The Combustion Institute, 1981
OF
TURBULENT
FIRE
PLUMES
FRANCESCO TAMANINI Factory Mutual Research Corporation, Norwood, Massachusetts 02062
The paper presents a model of fire plumes in which oxygen and fuel are allowed to be present in the average at a given height, even though chemical processes are assumed to be infinitely fast. The treatment of the composition fluctuations is derived from the k-~-g model of turbulence, in its adaptation to the case in which the variation of average properties across the plume is neglected (top-hat profiles). Predictions from the simple, one-dimensional integral model for an intermediate-scale propane flame show that important quantities, such as flame height and burning rate distribution in the flame, can be calculated w i t h accuracy close to that provided by the two-dimensional k-~-g approach. While final model validation still requires comparison with experimental measurements, which are not available at present, the model offers progress toward economical predictive techniques for simple combustion problems.
1. I n t r o d u c t i o n
Advanced models of turbulence and associated numerical techniques, ~ while justified for the purpose of improving our understanding of the details of turbulent flows, are probably not the most convenient tool for engineering calculations involving practical flows of simple geometry. Integral techniques have the potential for offering an alternate approach in which the important overall features of the flow are modeled realistically, at a significant saving in complexity of the formulation and in the time required to obtain solutions. In the integral approach such reduction in model complexity is achieved by assigning a prescribed shape to the cross-stream variation of flow properties; as a result, the model contains only ordinary differential equations, with longitudinal distance as the independent variable. Steward's work z on the modeling of turbulent fire plumes is an often referenced application of the integral technique to a situation in which both buoyancy and chemical reaction control flow development. The most serious limitation of Steward's formulation, as well as of others conceptually similar to it, 3-5 is the procedure used to predict the rate of turbulent combustion: instantaneous reaction with fuel of all 2'3 or of a constant fraction 4 of the entrained oxygen, or assumption of a stationary laminar flame front inside the turbulent plume. ~ The feature common to these three models of
turbulent combustion is the neglect of unmixedness. This is a phenomenon observed in turbulent reacting flows, where, due to the eddy-like structure of the flow, lumps of unburnt fuel appear interleaved with lumps of unreacted oxidizer. The overall rate of combustion is in this case determined by the reaction occurring at the periphery of the eddies and by the rate at which the total amount of available interface surface increases because of eddy breakdown caused by turbulence. This paper presents a proposal for improving the combustion-related aspects of the models in Refs 2-5: the proposal is an adaptation to the integral approach of the combustion model used in earlier work ~ with the k-~-g technique. The model will be discussed in its application to the case of axisymmetric fire plumes. In order to simplify the formulation, we will neglect the variations across the fire plume of time-averaged values of the flow properties (tophat profile assumption).
1081
2. M e a n F l o w E q u a t i o n s
2.1 Overall Continuity The rate of change of total mass flow, rh, with height, x, is given by: dm --
dx
= rn o
(1)
1082
T U R B U L E N T COMBUSTION
where rn o is the mass of ambient air entrained per unit height by the fire plume. Typically the entrained flow is calculated using an entrainment expression such as that proposed by Ricou and Spalding, ~ which reads: m' o - K(poW) '/2
(2)
where W is the vertical momentum of the fire plume at the height in question, and Po is ambient density. Equation (2) is conceptually the same as Morton's entrainment expression, 7 but for the added advantage that K is independent of the assumed velocity profile across the flow. Since we will refer to Eq. (2) later, we note that the entrainment coefficient K has been found to vary from .28 for forced jets ~ to .56 for thermal plumes. ~ Rather than relying on the approach represented by Eq. (2) we have chosen to calculate the entrained flow from: nio = ck~r(~ + ~,)
(3)
where I~ and p.~ are the laminar and turbulent components of viscosity and ck is a constant. The analytical solution of the laminar plume problem by Yih 9 indicates ck = 12. On the other hand, for the turbulent case the data by George et al, 8 when combined with the experimental estimate for p. made by Beuther et al, ~~ imply a value for ct close to 7. This latter value is the one used in our calculations. It should be noted that the fact that ck assumes a lower value in a turbulent flow than in the laminar counterpart has already been observed for the case of forced jets 'L'2 (ck = 5.5 and 8, respectively). As will be shown, the model allows one to estimate Ix~ and, through rt,, the intensity of the turbulent transport. Therefore, we feel that Eq. (3) can account in some fashion for the effect on entrainment due to incomplete development of the turbulence in the low part of the flame, while Eq. (2) obviously cannot. Another reason for our adoption of Eq. (3) in this work is the fact that the empirical coefficient, ck, seems to be less dependent than K on buoyancy effects: c k increases from 5.5 to 7, while K doubles in going from turbulent forced jets to thermal plumes. If Am o is the total entrained air flow and m, the initial mass flow, then the integrated form of Eq. (1) becomes: rh = m , + a m o .
(4)
2.2 Species Conservation The conservation of chemical species, i, can be written as:
d -(u dx
= rh', + u
(5)
where Y, and Y,.o are the average mass fractions in the plume and in the ambient, and m~ is the rate of production of species i per unit plume height due to chemical reaction. Species which do not participate in the chemical reaction, as in the case of nitrogen, have zero source (m', = 0); in the case of oxygen and fuel, however, m ' and th; are not zero for as long as unburnt fuel is present. Steward ~ assumed fixed stoichiometry and wrote (for Yr > 0): mo~ "' = srh; = -Yo, ohi'o
(6)
where u is the mass fraction of oxygen in the ambient and s is the oxygen/fuel stoichiometric ratio. Since initially the jet flow contains no oxygen (Yox.o = 0), Steward's assumption (Eq. (6)) implies that Yox - 0, for as long as Yr > 0 in the plume (see Eq. (5)). We will show in the following sections how we propose to allow for presence of oxygen and fuel at the same height in the plume, by accounting for the fluctuations in composition due to the turbulence. Let us now write Eq. (5) in a simpler form by introducing a Shvab-Zeldovich transformation: if oxygen and fuel combine in constant ratios, then the variable r = Yox - sYr
(7)
is a conserved quantity in that it is not created or depleted due to chemical reaction. The conservation equation for F is: d - - (rrh) = ron~;
(8)
dx
where F o = Yox.o is the ambient value of the composition parameter, F. By taking advantage of the definition of the total entrained flow used in Eq. (4), it is possible to integrate Eq. (8) for the case of uniform ambient conditions and write: F ni =
r,ni,
+ FoAm o
(9)
where F, = -sYs. ~ is the value of the composition parameter in the jet flow. For given initial and ambient conditions, Eqs. (9) and (4) can be combined to yield an analytical expression for F as a function of rh:
r = ro + ( r , - r o ) . , , / , a .
INTEGRAL MODEL OF TURBULENT FIRE PLUMES A similar expression applies to the case of the mass fraction of nitrogen (substitute Y,, Y,.~ and Y,.o for F, r~ and F o in Eq. (10)). Knowledge of the value of F is not sufficient to calculate Yo~and Ys; it is the task of the combustion model to determine the Yo~ and Ys associated with a given F. For the moment, we simply note that, once Yo~, Ys, and Y are known, the concentration of products can be obtained from: Y , . = 1 -Yo~ - Y e - Y..
(11)
since no other species need to be considered. 2.3 Energy Conservation The conservation of thermal energy requires that: d dx
(c.Tm) = q'~h+ c.Torho - q',
t12)
where T is the average temperature across the plume and cp is the specific heat (assumed constant) of the gases. The two quantities ti'~hand 0" are, respectively, the energy supplied by chemical reaction and that lost by radiative emission per unit plume height. The same procedure adopted in connection with the introduction of F, can be extended to the case of the enthalpy by defining a total enthalpy O as; O = c~T + QfYi
(13)
where Qf is the heat released by reaction of unit mass of fuel. The definition of Qf implies that: O'h = -Qsrh}
(14)
It is easy to show that the transport equation for O reads: d --
dx
(Ore) = Oom o -
(fr"
(15)
Equation (15) can be integrated to yield the following algebraic expression: Orb = O,rh~ + OoArho - &Jr
(16)
where Aqr is the total power radiated by the flame up to the height in question. Note that Eq. (16) would be identical to Eq. (9) if it wasn't for the presence of the radiation term, 5tj,: following the presentation of the combustion model, we will introduce a simple treatment of radiation which overcomes the complications associated with the difference between Eqs. (9) and (16).
1083
2.4 Momentum Equation The vertical momentum, W, of a plume in a coflowing stream of velocity U o, increases with height because of action of buoyancy: dW dx
= a,(po - p)~rb2 + ~hoUo
(17)
where a~ is the acceleration of gravity; p is the average density; and b is the radius of the fire plume. Note that the formulas introduced up to (but not including) Eq. (17) did not require a specification of the shape of the radial profiles of species, temperature, velocity, etc. inside the fire plume. A choice, however, is implicit in Eq. (17). While the model could accommodate other shapes, at the cost of an increase in complexity, we will restrict the discussion of the combustion model to the case of uniform (top-hat) profiles. In this case, the relationship for the plume radius is very simple: b = rh/X/'~pW
(lS)
3. Combustion M o d e l 3.1 The y-Plane So far we have been dealing with time-averaged quantities; having limited the scope of the analysis to such variables is the main cause for the lack of realism in Steward's model, z The behavior of turbulent flows cannot be understood unless appropriate consideration is given to the intrinsically unsteady, time-dependent nature of the flow. In other words, any turbulence model should account, in some more or less elaborate fashion, for the presence of fluctuations. Let us visualize a situation where, at a given height in the plume, the composition parameter fluctuates within a given range. The instantaneous values, ~, will be functions of time and position in the cross-section of the plume, and such that their time-space average is F. In order to evaluate the average values of species concentrations and density one needs two ingredients: 1) establishment of the relation between ~ and the corresponding instantaneous oxygen, fuel fractions (Yox, Ys), density (~) and temperature (/'); 2) determination of the probability of occurrence of all possible ~'s. If one assumes that chemical processes are infinitely fast, then fuel and oxygen cannot coexist (at the same time and at the same point); the relationship between ~ (or Yl) and z/ is easily obtained by remembering that by definition it is: r = Yox - s~s
(19)
1084
TURBULENT COMBUSTION
and that non-zero values of gox and #s correspond, respectively, to ~ greater and less than zero. With regard to the relationship between 5 and the instantaneous inert mass fraction, ~,, we note that the two normalized average variables (F F o ) / ( r , - ro) and (Y, - Y o)/(Y,.~ - Y,.o) obey the same differential equation and have identical boundary conditions. From this it follows that, for the two normalized variables, instantaneous values must also be identical to one another; the sought algebraic relation between ~ and #, is readily obtained. The concepts introduced here have already been used extensively in the literature: this is the reason for presenting them in cursory fashion. More details are available elsewhere. Determining the relationship between instantaneous temperature and ~ poses conceptual and practical problems, because of the presence of a term in Eq. (15), the radiated power q',, which has no counterpart in Eq. (8). If it is possible to assume that the radiated power is a constant fraction, • of the power produced by chemical reaction: q', = Xq'ch = - • Qcrh),
(20)
then one can define the total enthalpy as: O* = c T + (1 - •
(21)
and happily note that the term q'r no longer appears in the equation for O'. By invoking the same argument used in connection with ~,, it is then easy to obtain the algebraic relation linking t to Unfortunately, in the general case, where Eq. (20) may not be an accurate expression for t~, a unique relationship between t and ~, may not even exist. While a rigorous solution of this problem is beyond the capabilities of the approach described above, approximate solutions can be devised and have been used. la We are discussing this issue because we believe that it has implications in terms of the extent to which the modeling procedure can reliably accommodate more elaborate radiation treatments. Having made the point we now set it aside and assume, for the sake of the following application of the model, that Eq. (20) does indeed offer an accurate representation of the losses due to radiation. We will therefore calculate O ~ from Eq. (16) in which we have set Atlr = 0 and O = O ~ Once Ys is available from the statistics of the fluctuations, T can be obtained from Eq. (21). 3.2 Probability Distribution and Average
Quantities As stated previously, the evaluation of the mean values of concentrations and density requires
knowledge of the probability density function (PDF) of ~. The probability function H(~t) is defined in such a way that II (~) d~ is the probability of existence of concentrations in the range "~ to "~ + d~. We will follow an approach in which we assume that the PDF has a fixed shape; the width and the center of the distribution are determined from the_values of F and its mean square fluctuation, g = ~/2. The quantity g in turn is obtained from a transport equation as discussed in the following section. We choose the same form for the PDF already adopted in our earlier work, L'~ namely the shape given by two parabolae, which is reasonably close to the results from experiments and for which we have already worked out the necessary algebra. 14 The average value, q~, of a variable ~b(composition or density) can be evaluated through the integral:
i .,em~ "r
(22)
where 6(~) represents the analytical expression which relates instantaneous values of ili to those of ~,. In the general case, 6 (~) is a linear function, which is different for z/ greater or less than zero. Using Eq. (22) it is possible to calculate all average compositions, in particular Ys" Once Y~ is known, the burning rate per unit height is obtained from: d
m," = --dx (Y s,h)
(23)
3.3 Modeling of g-Equation Consider a mass flow of magnitude th and mean composition F (the composition parameter can be thought of as the mass fraction of an inert, since chemical reaction does not directly affect it). Let us further assume that the composition fluctuates i__~_ntime about its mean value, F, and define g = "t 2 as the magnitude of the mean square of the deviations from the mean. If one visualizes the flow as being an isolated stream, the only mechanism that can alter the value of g is the action of molecular diffusion processes, which will tend to reduce the intensity of the fluctuations until they approach zero. If, on the other hand, fresh fluid of uniform composition Fo is periodically added to the main stream, the level of fluctuations is going to increase for as long as F # F o. These two processes, molecular diffusion and entrainment of ambient fluid, constitute respectively the basic mechanisms of dissipation and production of g. We will examine now in some detail the practical implementation of these concepts by considering first the production of fluctuations associated with entrainment. An increment in mass flow from rn to rn + 8hi, due to entrainment of fluid of composition F o, brings
INTEGRAL MODEL OF TURBULENT FIRE PLUMES about an increase, 8g, of the level of fluctuations which can be calculated from:
1085
standard k-e-g model ~ as proposed by Spalding. ~" After some manipulation, the complete differential equation for g becomes:
rn (g + ~g) = [ g + (~r)~l -
-
d
E
. , ~-c~2 p -k- glrb 2 --(gin)= (F o - F ) 2too.
m + ~rh
(27)
dx
~m +(Fo-F-aF) ~
-
m+~m
(24)
The last term in the right-hand side of the equation represents the contribution to the overall level of fluctuations due to the fact that unit mass of flow contains the fraction ~rh/(rh + 8m), whose composition is Fo, while the mixture (after entrainment) has average composition F + ~F. The shift by 8F in the composition of the mixture, associated with entrainment of ~rh, is responsible for the second order increase (~F)2in the fluctuations of the original stream. The presence of mass ratios as weighing factors in Eq. (24) implies that g should be interpreted as a density-weighted quantity. Strictly speaking, all mean quantities introduced in the preceding sections should be thought of as Favre averages and not as Reynolds averages. However, a rigorous implementation of this fact would involve modeling of additional terms in the k- and ~- equations and, for this reason, we have not pursued it here (see 9 15 Bflger ). Equation (24) can be cast in differential form by dropping higher order terms; thus, the result for the case where all dissipation is neglected is:
dx
rho., = c~ ~r ~ ,
(28)
rather than the total entrained flow, as would be indicated by a more rigorous derivation of Eq. (27). ~ The distinction has some relevance only very low in the flame and we have made the change for reasons of convenience, to avoid numerical problems in the very early steps of the turbulent calculation. The evaluation of k, ~ and I~~is discussed next. 4. T u r b u l e n c e
Quantities
4.1 Turbulence Kinetic Energy and Its Dissipation The transport equation for the turbulence kinetic energy is given by:
d(krh) - -
d --
where k and e are the turbulence kinetic energy and its dissipation; while c,2 is a constant. Note that we have introduced in the source term the flow entrainment associated with the turbulent fluctuations:
= Pk + G~ - ~:~
(29)
dx (g,n) = ( r o - r ) ~ , n o .
(25)
In the absence of molecular diffusion, Eq. (25) is an exact mathematical representation of the mechanism governing the growth of g. It is worth noting that, in the case of a jet of uniform composition, F~, issuing in an environment of uniform composition, Fo, Eqs. (8) and (25) can be solved analytically and the solution is: g = (ro - r)(r - r,).
(26)
The level of fluctuations given by this solution is the maximum attainable and corresponds in the PDF plane to two delta functions at ~/ = Fo and zt = F, and no intermediate compositions. As has been pointed out, the rate of production of g can be calculated from an expression (Eq. (25)) which is exact. For the case of the dissipation of g, however, the situation is not as fortunate and we must resort to modeling. While recognizing that this point may be the subject of future improvements, we will for the moment simply extend to the plug flow considered here the formulation used in the
where Pk and G~ are total production of k due to shear and buoyancy effects respectively, and ~ is total dissipation. By analogy with the exact expression obtained for the term representing production of g (cf. Eq. (25)), we write: 1
Pk = - - (U - U o ) 2 mo., .,
2
(30)
where U is the average vertical velocity of the fire plume and U o is the vertical velocity of the free stream (in the case considered here it is Uo = 0). For the production term representing generation of k due to buoyancy we write: 17
Gk = agco X / ~ k ~rb2
(31)
*M. A. Delichatsios has brought to our attention the fact that Eq. (27) can be obtained by direct manipulation of the two-dimensional boundarylayer form of the model equation for g. His contribution is gratefully acknowledged.
T U R B U L E N T COMBUSTION
1086
where c o is a constant (=.15) and ~ is the root mean square value of the density fluctuations. Finally, the total dissipation is given by: E~ = p~Trb ~
(32)
For the dissipation ~, we follow standard practice L~8-2oand write:
- dx
-
{c,~P k + c,~G k - c,2Ek}
(33)
k
where c~, c~3 and c,2 are constants. We will later see that values for these constants, close to those considered to be "standard" L~8-2onamely c,~ = 1.44, c,2 = 1.82 and c,3 = .95, lead to reasonable model predictions.
TABLE I Properties of propane flame used as test case Jet velocity, u~ (m/s) Nozzle radius, b~ (m) Fuel jet temperature, T t (K) Ambient temperature, T O (K) Specific heat, cp (J/kg/K) Oxygen/fuel ratio, s (-) Heat of combustion, Qr (J/kg) Radiated fraction, • (-) Ambient oxygen mass fraction, Yo~.o (-) Jet fuel mass fraction, Yr (-)
3.251 .00635 300. 300. 1200. 3.636 4.64 ::< 1 0 7 .26 .233 1.
5. Model Predictions 4.2 T u r b u l e n t Viscosity The entrainment formula given in Eq. (3) contains the turbulent viscosity, p,. We obtain this quantity from the known expression, used by the k-~-g model, ts which reads: k2
p~ = c~ pV/-~op - -
(34)
s
where c = .09 is a constant. Note that the average density ~pop)~/2 is used instead of p. At this point the model is complete and its application to a specific situation requires only the choice of appropriate initial conditions. 4.3 Initial Conditions As we have indicated in previous applications of the k-e-g technique, L~3 the problem of how to start a turbulent flow calculation has yet to find a satisfactory solution. The reason is that the models used apply to fully developed turbulence and cannot deal with the initial growth in the transition region. However, the arbitrary nature of the procedures usually adopted to introduce turbulence in the flow should not justify incomplete documentation. In the calculations presented here, we have assumed initial values for k, e and g right at the exit of the jet nozzle as follows: k = 10-au~ g=0
(35) (36)
and for ~ the value that corresponded to equilibrium (Pk = Ek, cf. Eqs. (30) and (32)). The transport equations for k, e and g (Eqs. (29), (33) and (27)) then dictated the subsequent development of these quantities.
5.1 Constants Selection The model used in this work contains a number of adjustable constants, whose values can be chosen to improve agreement between theory and experiment. Of special concern to fire research is the variation with height in the plume of the volumetric burning rate. This quantity is important because convective and radiative heat release and, therefore, flow development are all linked to the location and the intensity of the chemical reaction. Since no direct measurements are available on the burning rate distribution in buoyancy-controlled flames, we have used as a point of reference the numerical predictions for the largest of three intermediate-scale propane flames, considered in our earlier work with the k-e-g technique. ~ The properties of the propane jet flame are listed in Table I and the model constants finally adopted are given in Table II. While the values of constants c,~(=1.44), c,2 (=1.82) and c.(=.09) are the same as those used in two- or three-dimensional applications of the k-e-g approach, changes in the other constants were made necessary by the particular nature of the integral approach as implemented here. So, for example, in
TABLE II Model constants Entrainment formula (Eq. (3)): g-equation (Eq. (27)): k-equation (Eq. (31)): e-equation (Eq. (33)):
Turbulent viscosity formula (Eq. (34)):
c k = 7.
ca2 = cp = c,1 = c,2 = cE~ =
.8 .15 1.44 1.82 .95
c~ = .09
INTEGRAL MODEL OF T U R B U L E N T FIRE PLUMES order to obtain a flame height in agreement with that given by our earlier predictions,~ (which in turn was found to agree with experiment), we had to reduce the value adopted for %2 down to .8. The values for c 2 quoted in the literature are in the 19 g 12 20 range 1.25 to 2.0 to 2.8 . For our earher flame predictions ~ we set c~z = 1.82. To counteract the effect of the resulting high level of density fluctuations on G,, the buoyancy production of turbulence kinetic energy, we have set c = .15. This choice should be seen against the fact that the true value of the correlation coefficient c , based on thermal plume dataY is close to .6-.7. The choice of c,3 = .95 was dictated by the desire to obtain above the flame tip the same entrainment coefficient (K = .56) indicated by thermal plume data. s
12
,
1087 --
L o l ~
9 t,,,.. 1/*? y\
1500
13oo
,,oo
Height, x(m) Fie. 2. Predictions of the integral model of average fuel and oxygen mass fractions, and temperature at various heights in the p r o p a n e flame.
5.2 Numerical Solution The solution of the model equations is quite straightforward. There are five ordinary differential equations (Eqs. (1), (17), (27), (29) and (33)) for rh, W, g, k and e, and four algebraic expressions (Eqs. (3), (10), (16) and (34)) for rh~,, F, O* and Additional relationships are used for the statistical evaluation of mean flow properties. The flame calculation described in the next section required less than 150 integration steps in the x-coordinate and less than 4 sec of CPU time on IBM 370/158.
Is~.
5.3 Flame Calculation A comparison between the predictions of variation of burning rate with height given by the full k-e-g
.ool2 -r c
,' ,;
9 ',,
Propane flame UI= 5.251 rn/s
0008
g:.ooo,~ m
,4
.8 Height, x(m)
1.2
1,6
Flc. 1. Calculated variation with height of the burning rate per unit plume height. The result of a calculation with the integral model is compared with the prediction from the standard k-(-g approach.
model ' and by its adaptation to the integral approach are shown in Fig. 1. As can be seen, the trends followed by the more complex formulation are faithfully reproduced by the simpler model developed in the present work. Even better agreement could have been obtained by additional model tuning (read: improved selection of constants). This further refinement didn't appear to be worth the trouble, since final validation of the model can only come from experimental data. Figure 2 shows the variation of oxygen and fuel mass fraction, temperature and unburnt fraction of jet fuel flow as a function of height. As can be seen, very near the base of the flame oxygen starts to appear in the fire plume, long before fuel is consumed. Furthermore, the average temperature of the fire plume reaches a maximum value of about 1450 K, which is well below the adiabatic flame temperature for this system (2023 K). The ratio between the flux of unburnt fuel (mYe) and its value at the nozzle (rh~) varies gradually from 1, at x = 0, to 0 at the flame tip. As can be seen in Fig. 2, 99% of the fuel is consumed at the height of 1.16 m. We are currently developing a novel experimental technique which will enable us to make a direct measurement of this quantity. Calculated values for vertical velocity, U, plume radius, b, and entrainment coefficient, K, are shown in Fig. 3 as a function of height. The vertical velocity first drops quite rapidly to about half the nozzle value, then decreases more gradually in the remaining part of the flame. The plume width increases steadily and so does the entrainment coefficient, which tends to approach asymptotically the thermal plume value (K = .56). The curve for K in Fig. 3 and that for the burning rate in Fig. 1 show an initial rapid rise, which is then followed by an inflection point. This occurs in the low part of the flame, where the turbulence undergoes rapid growth
1088
TURBULENT COMBUSTION
bU .2
K
Nomenclature
4
.155~N~
ag
Ent~
"4
b Cg2 r
acceleration of gravity, (m/s 2) fire plume radius, (m) constant in g-equation (Eq. (27)), (-) constant in entrainment formula (Eq.
(3)), (-) Cp C~l, C~2, C~3 C~ r
.4
Hegixh(m t, ) .8
1.2
i
g Gk
FIG. 3. Calculated variation with height of mean
vertical velocity, plume radius and entrainment coefficient as given by the integral model. in a region of numerical transition. We have verified to our satisfaction that this anomaly is not caused by numerical problems with the solution, but is dictated by the model equations.
k K m
m; Pk
q'~ 6. C o n c l u s i o n s
We have presented an integral model of fire plumes which uses the approach to the modeling of turbulent processes introduced by the k-e-g technique. Even though on the conceptual front our adaptation of the k-~-g technique offers few new ideas, the importance of this work lies in the fact that the integral model is shown to yield results, similar to those obtainable from the full 2-dimensional k-~-g approach, but at a fraction of the computer time (4 sec instead of 150 sec of CPU time on IBM 370/158 for. a flame calculation). This reduction in model complexity, with little sacrifice in terms of accuracy, is necessary for applications which require several flame calculations such as, for example, in problems of fire spread. Validation of the model will require comparison with experiment; as already mentioned, we are working on a facility in which we will be able to measure the burning rate distribution in the fire plume, as well as other important parameters. We can already anticipate that several aspects of the model will need improvements; the treatment of entrainment is one of them, particularly if one intends to use this formulation to calculate pool fires. The advantage of this approach lies in its simplicity, as demonstrated by the relative shortness of the computer program which implements it. The same computer program can easily be modified by external users to incorporate model expressions different from those described here.
q', Q~
specific heat, (J / kg/K) constants in e-equation (Eq. (33)), (-) constant in turbulent viscosity formula (Eq. (34)), (-) constant in formula for buoyancy-generated turbulence (Eq. (31)), (-) mean square fluctuations of composition parameter, ( - ) production of k due to buoyancy (Eq. (31)), (kg m / s 3) turbulence kinetic energy, (m2/s ~) entrainment coefficient (Eq. (2)), (-) total mass flow, (kg/s) rate of production of species i per unit height, (kg/m/s) ambient fluid entrained per unit height, (kg/m/s) production of k due to shear (Eq. (30)),
(kg m/s 3) rate of release of chemical energy per unit height, (W/m) radiated power per unit height, (W / m) heat of combustion of unit mass of fuel,
(J/kg) $
t,T U W X
y,,Y, "y,F hrho e
oxygen ] fuel stoichiometric ratio, (-) temperature, (K) mean vertical velocity, (m/s) vertical momentum, (kg m / s ~) longitudinal distance, (m) mass fraction of species i, (-) composition parameter (Eqs. (7) and (19)), (-) entrained mass flow, (kg/s) dissipation of turbulence kinetic energy, (m2/s 3) total dissipation of k (Eq. (32)), (kg
m/s ~) 0 09
total enthalpy (Eq. (13)), (J/kg) corrected total enthalpy (Eq. (21)),
p. P.T II p •
laminar viscosity, (kg/m/s) turbulent viscosity, (kg/m/s) probability density function of ~, (-) density, (kg/m 3) radiated fraction, (-)
(]/kg)
Superscripts -
~ '
indicates average value indicates instantaneous quantity indicates fluctuating component
INTEGRAL M O D E L OF TURBULENT FIRE P L U M E S
Subscripts f n ox pr 0 1
fuel diluent (nitrogen) oxygen products ambient conditions fuel jet conditions
Acknowledgments The author is grateful to John de Ris for useful discussions. The work was carried out under joint sponsorship from Factory Mutual Research Corporation and the National Bureau of Standards under Grant NBS G7-9011,
REFERENCES 1. TAMANINI,F., Combustion and Flame, 30, 85-101 (1977). 2. Steward, F., Combustion Science and Technology, 2, 203-212 (1970). 3. NIELSEN,H. J. ANDTAO, L. N., Tenth Symposium (International) on Combustion, pp. 965-972, The Combustion Institute, 1965. 4. BRZUSTOWSKI,T. A., in "Turbulent Combustion," Progress in Astronautics and Aeronautics, Vol. 58, L. A. Kennedy editor, pp. 407-430, American Institute of Aeronautics and Astronautics, 1978. 5. WILCOX, D. C., AIAA Journal, 13, 3, 381-386 (1975). 6. RIcou, F. P. AND SPALDING, D. B., J. Fluid Mechanics, 11, 21-32 (1961). 7. MORTON,B. R., Tenth Symposium (International) on Combustion, pp. 973-982, The Combustion Institute, 1965. 8. GEORGE, W. K., ALPERT, R. L. AND TAMANINI, F., Int. J. Heat Mass Transfer, 20, 1145-1154 (1977).
1089
9. YIH, C. S., Proceedings of the 1st U.S. National Congress of Applied Mechanics, E. Sternberg editor, pp. 941-947, American Society of Mechanical Engineers, New York, 1951. 10. BEUTHER, P. D., CAPP, S. P. AND GEORGE, W. K., "Momentum and Temperature Balance Measurements in an Axisymmetric Turbulent Plume," ASME paper 79-HT-42 presented at the Joint A S M E / A I C h E 18th National Heat Transfer Conference, San Diego, California, August 6-8, 1979. 11. SCHLmHTJNG,H., "Boundary Layer Theory," p. 184 and p. 609, Mcgraw-Hill, 1960. 12. SPALDING, D. B., "'Combustion and Mass Transfer," p. 207, Pergamon Press, 1979. 13. TAMAN1NI,F., Seventeenth Symposium (International) on Combustion, pp. 1075-1085, The Combustion Institute, 1979. 14. TAMANINI,F., "Numerical Model for the Prediction of Buoyancy-Controlled Turbulent Diffusion Flames (A Description of the F U E L J E T 2 Computer Program)," Factory Mutual Research Technical Report, Serial No. 22360-1, RC-B-47, October 1975. 15. BILGE8, R. W., Prog. Energy Combust. Sci., 1, 87-109 (1976). 16. SPALDING,D. B., Chemical Engineering Science, 26, 95-107 (1971). 17. TAMANINI,F., "On the Numerical Prediction of Turbulent Diffusion Flames," paper presented at the Joint Meeting of the Central and Western States Section of the Combustion Institute, Southwest Research Institute, San Antonio, Texas, April 21-22, 1975. 18. LAUNDER, B. E. AND SPALOINr D. B., "Mathematical Models of Turbulence," Academic Press, 1972. 19. GIBSON, M. M. ANt) LAUNt)Ea, B. E., ASME ]. of Heat Transfer, 98, 81-87, 1976. 20. TAMANINI, F., ASME J. of Heat Transfer, 100, 659-664, 1978.
COMMENTS S. Galant, Societe Bertin, France. In largely reducing computational times, your work allows for easy sensitivity analysis of the results of the model parameters. Have you carried out such an analysis for the k-~-g model parameters? If so, how do they affect numerical results? Author's Reply. We have not carried out a systematic analysis of the effect of changes in the values of the various constants. The discussion given in the paper on the criteria followed in selecting such
values, indicates that we have fine-tuned the model only to a limited extent. The flame height depends rather strongly on the value of ce~ and we have taken advantage of this fact to fit the integral model to the results from the 2D calculations. Any additional model optimization will have to be done by comparison with experimental data. I should point out that, in any case, this process of optimization cannot use all the degrees of freedom given by the seven (7) model constants. In fact, it would be hard to justify the introduction of drastic changes, which
1090
TURBULENT COMBUSTION
w o u l d m a k e the set o f c o n s t a n t s m u c h d i f f e r e n t from the o n e generally a c c e p t e d a m o n g k-c-g modelers.
fire p l u m e are very weakly d e p e n d e n t on the s h a p e of the P D F . T h i s c o n c l u s i o n h o l d s also for the integral m o d e l p r e s e n t e d in the paper.
H. Milter, Harvard University, USA. D i d y o u assume a s h a p e for y o u r P D F , or solve a P D F e q u a t i o n ? If the former, t h e n of c o u r s e y o u have in fact more t h a n the 7 p a r a m e t e r s y o u appear to have; y o u r results will p r e s u m a b l y d e p e n d on the a s s u m e d P D F . C a n y o u say, at this time, w h a t is the s e n s i t i v i t y of the results to the a s s u m e d s h a p e ?
Author's Reply. As i n d i c a t e d in the paper, we a s s u m e that the P D F in the 7y-plane is u n i q u e l y d e f i n e d by each pair of v a l u e s for F a n d g. O u r p r e v i o u s work w i t h t h e 2I) v e r s i o n of t h e k-e-g m o d e l (Ref. 1) h a s s h o w n that the overall features of the
1. Creighton, Lawrence Livermore Laboratories, USA. F r o m t h e oral p r e s e n t a t i o n , it w a s n o t clear to m e w h i c h c o m p u t e r code or n u m e r i c a l m e t h o d was u s e d to calculate the K-c-g model.
Author's Reply. T h e details of o u r 2 D calculations u s i n g the k-~-g m o d e l are reported in Ref. 1. T h e c o m p u t e r code a d o p t s the s o l u t i o n a l g o r i t h m for parabolic flaws d e v e l o p e d b y S p a l d i n g a n d coworkers at I m p e r i a l College a n d k n o w n u n d e r the code n a m e of G E N M I X .