SANDIA REPORT

N U REGICR -4855

SAN D8 7-8 203

R3

Printed July 1 987

Development and Application of a

Computer Model for Large-Scale

Flame Acceleration Experiments

K. D. Marx

Prepared by

Sandia National Laboratories

Albuquerque, New Mexico 718and Livermore, California 94550

for the United States Department -ofiEnergy

under Contract DE-AC04-76Dbpb0089

ilk

A A

'I

II1l~ VI

'iL>iI !H TI _

It'

SF2900Q(8-811

Y11

I'
L

Issued by Sandia National Laboratories, operated for the United States

Department of Energy by Sandia Corporation.

NTICE: This report was prepared as an account of work sponsored by

an agency of the United States Government. Neither the United States

Government nor any agency thereof, nor any of their employees, nor any

of the contractors, subcontractors, or their employees, makes any warranty, express or Impiied, or assumes any iegai liability or responsibility

for the accuracy, completeness, or usefuiness of any Information,

apparatus,, product, or process disclosed, or represents that its use would not Infringe privately owned rights. Reference herein to any

specific commercial product, process, 'or service by trade name,

trademark, manufacturer, or otherwise, does not necessarly constitute

or Impiy its endorsement, recommendation, or favoring by the United

States Government, any agency thereof or any of their contractors or

subconractors. The views and opinions expressed herein do not

necessarly state or refiect those of the United States Government, any

,agency thereof or any of their contractors or subcontractors.

Available from

Superintendent of Documents

U.S. Government Printing Office

Post Office Box 37082

Washington. D.C. 20013-7982

* and

National Technical Information Service

Springfield, VA 22161
NUREG/CR-4855

SAN D87-8203

R3

DEVELOPMENT AND APPLICATION OF

A COMPUTER MODEL FOR LARGE-SCALE

FLAME ACCELERATION EXPERIMENTS

K. D.Marx

July 1987

Sandia National Laboratories

Livermore, CA 94550

Operated by

Sandia Corporation

for the

U. S. Department of Energy

Prepared for

Division of Reactor Accident Analysis

Office of Nuclear Regulatory Research

U. S. Nuclear Regulatory Commission

Washington, DC 20555

Under Memorandum of Understanding DOE 40-550-75

NRC FIN No. M1246
A
DEVELOPMENT AND APPLICATION OF

A COMPUTER MODEL FOR LARGE-SCALE

FLAME ACCELERATION EXPERIMENTS*.

K. D. Marx

Thermofluids Division

Sandia National Laboratories

Livermore, CA 94550

ABSTRACT

A new computational model for large-scale premixed flames is developed and applied to

the simulation. of flame acceleration experiments. The primary objective is to circumvent the

necessity for resolving turbulent flame fronts; this is imperative because of the relatively coarse

computational grids which must be used in engineering calculations. The essence of the model

is to artificially thicken the flame by increasing the appropriate diffusivities and decreasing the

combustion rate, but to do this in such a way that the burn velocity varies with pressure, temperature, and turbulence intensity according to prespecified phenomenological characteristics. The

model is particularly aimed at implementation in computer codes which simulate compressible

flows. To this end, it is applied to the two-dimensional simulation of hydrog .en-air flame acceleration experiments in which the flame speeds and gas flow velocities attain or exceed the speed

of sound in the gas. It is shown that many of the features of the flame trajectories and pressure

histories in the experiments are simulated quite well by the model. Using the comparison of experimental and computational results as a guide, some insight is developed into the processes

which occur in such experiments.

*This work was performed at the Combustion Research Facility and supported by the

U. S. Nuclear Regulatory Commission.

ii i/iv

Table of Contents

Acknowledgements.... ...... ...... ...... ...... ...... ...... ...... xv

Executive Summary .. .. ...... ...... .. .. ...... ...... ... ... ... .... XV

I. Introduction .. .. .... ...... .......... .. .. .... ...... .......... 1

III. Computer Code and Turbulence Mode! .. .. .. ...... ...... ...... ..... 5

Flow Variables and Chemistry.. . ..... .. ...... ...... .......... 5

The k-c Turbulence Model .. .. .... ...... ...... ...... ........ 7

Ill. The Combustion Model. .. .. ...... .... ............ .. .... ...... 10

Properties of Turbulent Flames and Numerical Flames. .. .... ........ 10

Development of a Model With Appropriate Properties .. .. .... ........ 13

Dependence on Pressure and Preflame Temperature.. .. .. .. ........ 15

Summary of the Model .. .. ...... ...... .. I.... ...... ......... 17

IV. Comparison of Computational Results With Experiment .. . ..... .. ........ 20

Experimental Considerations .. .. .. ...... ...... ...... ........ 20

Comparison of Flame Trajectories and Pressure Histories .. .. .. ...... 22

Further Discussion of the Computational Results. ..... .. .. .. ...... 35

Results of Variations in the Computational Parameters. .. .... ........ 40

V. Concluding Remarks. .. .. ...... ...... ...... ...... ...... ....... 45

Appendix A. Integration of Prod uction- Dissipation Terms .. . .. ....... 51

Appendix B. Difference Formula for c at Walls. .. .. ...... ...... ...... .... 54

References .. .. ...... ...... ...... ...... ...... .... ...... ........ 55

V/Vi

List of Illustrations

.1 Sketch of the FLAME Facility .................

2 Schematic of the FLAME Facility with obstacles installed. ...

3 Problem domain and finite-difference grid. .. .. ...... .. ....

4 Length scales in the flame .. .. .... ........ ...... .. ....

5 Turbulent burn velocity vs turbulence intensity.... .. .. .. ....

6 Arrival time data for Experiment F-22.... .. .. .. ........

7 Flame front profiles for Experiment F-22 .. .. ......... .. ....

8 Arrival time data for Experiment F-23.... .. .. .. .......

9 Flame front profiles for Experiment F-23.... .. .. .. .. ......

10 Computed and experimental flame trajectories (Experiment F-23)

11 Computed and experimental pressure histories at 14.6 m

(Experiment F-23, transducer PKU14).. .. .... ....... .. .. ..

12 Computed and experimental pressure histories at 18.6 m

(Experiment F-23, transducers PKU5 and PKU1O)... .. .. ...... 13 Computed and experimental pressure histories at 21.5m

(Experiment F-23, transducers PKU6 and PKU9)... .. .. .. ....

14 Computed and experimental pressure histories at 25.6 m

(Experiment F-23, transducers PKU7 and PKU18)... .. .. .. ....

Paqe

2

3

8

. .. . 16

. .. . 23

. .. . 24

. .. . 25

26

. .. . 29

15 Three-dimensional plots of computed pressure (baseline parameters). ...

16 Three-dimensional plots of computed temperature (baseline parameters) . ...

17 Three-dimensional plots of computed turbulent kinetic energy

.I(baseline parameters)............. . .. . .. .... .. .. .. .. .. .. .. .

18 Three-dimensional plots of computed combustion rate

(baseline parameters). .. ...... ............. ....... . . .. .. ....

19 Flow pattern in two chambers at t =310 ms....... . .. .. .. .. .. .. ....

20 .Flame trajectories, including computation on fine grid ............

21 Pressure histories, including computation on fine grid ............

22 Three-dimensional plots of pressure computed on fine grid.... .. .. .. ....

23 Flame trajectories, including computation with altered model... .. .. .. ....

24 Pressure histories, including computation with altered model... .. .. .. .... 25 Pressure histories, including computation with n = 2.... .. .. .. .. ......

30

31

32

33

36

37

38

39

41

42

43

44

46

47

48

Vii/viii

List of Tables

Page

1. Results of Tests of the Combustion Model. .. .... ...... ...... ........ 19

I11. Summary of FLAME Tests With Obstacles* . . . . . . . .. .. . . . . . .. . . . . . 21

'111. Summary of Experimental and Numerical Parameters .. .. ......... .. .... 27

IV. Location of Pressure Transducers in Experiment F-23. .. .... ...... ...... 34

ixlX

List of Symbols

Roman:

A

B

cp

CPO

C1

c2

CM

D

E (App. A only)

f (T)

fp(q,k,c) (App. A only)

h(p, T.)

I

k

km

kT

L

Le

Lm

m

n

p

P (App. A only)

Pr

P' (App. A only)

Q (App. A only)

r

RT

Sc

5 Ti

ST (Table I only)

Coefficient in burn velocity formula

Combustion coefficient

Coefficient in burn velocity formula

Specific heat at constant pressure

Constant in approximation for specific heat.

Slope in approximation for specific heat

Progress variable

coefficient 'in production of c

Coefficient in dissipation of E

Coefficient in eddy viscosity formula

Diffusion coefficient

Special burn rate function

a;: Vu/g

Special burn rate function

Spe cific Enthalpy of Z~th species

Specific energy

Turbulent kinetic energy

Effective turbulent kinetic energy

Threshold turbulent kinetic energy

Thermal conductivity

Turbulent length scale

Effective turbulent length scale

Threshold turbulent length scale

Pressure exponent in burn velocity formula

Molecular weight .of ith species

Temperature exponent in burn velocity formula

Pressure

or: VU

Prandtl number

Clufp

pk

Rate of energy release due to chemistry

Thermal expansion ratio

Turbulent Reynolds number

Schmidt number

Laminar burn velocity

Turbulent burn velocity

Initial turbulent burn velocity

Expected turbulent burn velocity

dST /dVlik

xi
t

tt

tT

T,

tx

U

U,.

U1

Vf 0

VT

w

y (App. B only)

Yi

z

Time

Typical transit time

Turbulence time scale

Temperature

Adiabatic flame temperature

Initial gas, temperature

Unburned gas temperature

Thermal diffusion time

.Gas velocity

Wall shear velocity

Turbulence intensity

Initial flame speed

Turbulent flame speed

Channel width

Spanwise coordinate

Vertical coordinate

Coordinate measuring distance from wall,

Mass fraction of the ith species

Mass fraction of product species for completed burn

Lengthwise coordinate

Greek:

6L

At

AX

Ay (App. B only)

ALT

Vi

ý (App. A only)

p

Ai

Refers to 3-transformation

Unit tensor

Computational flame thickness

Expected computational flame thickness

Laminar flame thickness

Turbulent flame thickness

Time step

Grid spacing

Grid spacing

Turbulence dissipation rate

Von K~rm~n constant

Total viscosity

Laminar viscosity

Eddy viscosity

Stoichiometric coefficient of jth species

OlE

Total density,

Density of it" species

Xii
or

1*

Rate of consumption of fuel

Rate of change of density of it species due to reactions

Stress tensor

Coefficient in diffusion of E

Wall shear stress

Thermal diffusivity

Subscript :

e

f

i (App. B only)

L

T

Effective values

Fuel

Species

Grid index

Laminar

Product species

Turbulent

Refers to turbulent dissipation rate

Superscripts:

base (Table I only)

C

e (Table I only)

T

prime (') (on ci,, ST)

prime (') (on A~, tT, AT)

Refers to baseline calculation

Chemistry

Expected values

Transpose of tensor

Differentiation

Values of quantities in final form of combustion model

Abbreviations:

CJ

DDT

FLAME

Chapman-Jouguet

Deflagration-to-detonation transition

Flame Acceleration Measurements and Experiments facility

xiii/xiv

Acknowledgements

Many thanks are due M. P. Sherman, S. R. Tieszen, J. T. Hitchcock, J. E. Shepherd, J.

H. S. Lee, W. B. Benedick, P. K. Barr, and A. R. Kerstein for valuable discussions of this project.

I would also like to thank P. Worthingt on, C. Tinkler, M. Berman, and B. R. Sanders for. their

continued advice and support.

xv/xvi
i
Executive Summary

This report represents a major portion of the work performed at Sandia National Laboratories, Livermore under the Sandia Hydrogen Program for the U.S. Nuclear Regulatory Commission. Computational methods are described which permit techniques previously used for

modeling small-scale laboratory experiments to be extended to the simulation of experiments

with much larger length scales. The specific configurations considered are the propagation of

hydrogen-air burns in the FLAME facility at Sandia National Laboratories, Albuquerque. The

largest length scales appearing in these experiments approach in size those which might occur

in nuclear reactor accident scenarios.

These length scales. exceed the thicknesses of the hydrogen flames by a large factor.

Because of this, it was necessary to invent a new procedure to represent the flow variables on

a finite-difference grid in which computer limitations require that the grid spacing exceed the

physical flame thickness. The new method consists of imposing minimum values of turbulence

intensity and turbulent length scale so as to artificially thicken the flame. This is done in such a

way that the burn velocity obeys established dependencies on turbulence level, pressure, and

unburned gas temperature.

The resulting model is applied to the simulation of one of the experiments performed with

obstructions in the channel of the FLAME facility. By using only the initial flame velocity to adjust

parameters in the model, a quite successful representation of the experimental flame trajectory

and pressure histories is obtained. The calculation does display some inadequacies; it fails to

reproduce a sudden acceleration midway in the experiment, and it exhibits a detonation near

the end. Nonetheless, it provides a meaningful study of many features of the gas flow and flame

propagation.

As one example of this, the computations point out the way in which the acceleration of

the flame and its propagation in the choking regime are accompanied by large increases in the

amplitude of pressure waves which reflect off the obstacles. This mechanism contributes to the

positive feedback which drives the flame acceleration in two ways: (1) it leads to an increased

combustion rate via compression and shock heating and (2) it produces increased turbulence

through the generation of large flow velocities with high shear, thereby augmenting the thermal

mixing. These increasing pressure and turbulence levels are also undoubtedly indicative of a

buildup of sensitivity to detonation. This emphasizes the importance of a careful interpretation

of experiments with regard to the degree of confinement of the gas, which has been wellestablished experimentally. The results are consistent with a reduction in flame acceleration in

the presence of pressure relief due to venting. (it is noted, however, that there is a potential

for increased production of turbulence in the presence of partial venting which remains to be

investigated.)

In summary, this work contains two major contributions: (1) the development of the numerical model, and (2) the use of the computer calculations to provide insight into the processes

accompanying flame acceleration in experiments in the FLAME facility. It defines the attributes

and limitations of current computational capabilities in this research area, and suggests additional steps that could be taken to improve them.

xvii/xviii

I. INTRODUCTION

In this paper, we describe a computational model for simulating the propagation of premixed turbulent flames which extend over spatial length scales which are large compared to

those which characterize the turbulent flame zone. In these situations, the, geometric complexity

is often such that it is impossible to provide accurate resolution of the flame front on a computational grid Which encompasses the entire domain of interest. The specific configurations to

which this work has been applied are experiments performed in the FLAME facility' - 3 at Sandia

National Laboratories in Albuquerque, New Mexico (see Figure 1). The FLAME facility consists

of a reinforced concrete channel 2.44 m high by 1.83 m wide by 30.5 mn long (8 ft by 6 ft by 100

ft). The particular experiments considered here involve the ignition of a hydrogen-air mixture at

one (closed) end and the propagation of the flame past periodically spaced obstacles toward

the other (open) end (see Figure 2).

Turbulence- generated ahead of the flame as the unburned gas flows past the obstacles

results in an increase in the combustion rate. This, in turn, causes the flow velocity of the

unburned gas to increase, with a concurrent increase in the turbulence level. Hence, there is a

positive feedback mechanism which leads to acceleration of the flame as it propagates down

the channel.

.This mechanism for flame acceleration has been described previously for experiments on

a smaller scale.4 -

6 Some work has also been done on large-scale experiments .7 The primary

purpose of the FLAME facility is to provide relatively large-scale results for hydrogen-air burns.

The ultimate goal of the research is to determine scaling laws which will permit assessment

.of potential hazards to nuclear reactors in case of hydrogen production in reactor accidents

(see,. e.g., Reference 8). The length scales existing in the FLAME facility are much larger than

those typically occurring in the small-scale experiments described in References 4-6, but. are

somewhat smaller than some which might be encountered in a reactor containment.

The purposes of this paper are to describe our method for simulating flame acceleration in such large-scale configurations, to compare computational results with a limited set of

experimental data, and to provide some analysis of the processes occurring in the experiments.

The computer modeling of reacting flows as complex as those existing in the present

experiments is extremely difficult. There are, in principle, a large number of chemical reactions

to be accounted for. In all but the simplest flow configurations, one must approximate the true

chemical kinetics with a reduced model. Furthermore, as noted above, turbulence plays a crucial

role in determining the rate of mixing of burned and unburned gas, and therefore has a great .influence over the rate of combustion. A complete solution to the Navier-Stokes equations for

turb 'ulent flow is manifestly impossible in practical engineering calculations such as these. It is

therefore necessary to provide a turbulence model which permits the computation of turbulence

intensities and time scales. In this work, we have used the k-E model. 9-10 Given such a method

for the determination of turbulence parameters, one must couple it to the combustion model in

order to have a realistic determination of the rate at which the unburned gas is consumed. There

have been various proposals for reasonable waysito do this. Most notable for our purposes are

the Magnussen-Hjertager model'" and the eddy breakup model."- 1

3 .

1
Closed end of channel

Figure 1. Sketch of the FLAME facility. (From Reference 3.) The vertical coordinate is assumed ignorable in the two-dimensional computer calculations. (Note: FLAME is an acronym

for FILame Acceleration Measurements and Experiments.")

2
/STEEL LATTICE FRAME /PLYWOOD END DOOR /PLYWOOD OBSTACLE

REINFORCED CONCRETE

/VIEWPORT THERMOCOUPLE

RAKE

I~o~s~m I1.83m-*

7

0.74m t

1.83m

1i

NO

30.48m

STEEL REINFORCEMENT

PLATE I

2.44m

IFigure 2. Schematic of the FLAME facility with obstacles installed. (From Reference 34.) The

obstacles used in the experiments described in this work are shown to scale; the blockage ratio is !. The enclosed area in the top view is the two-dimensional domain of the computer 3.

calculations.

3
The Magnussen-Hjertager model is useful in situations where the flame propagation is

determined by the rate of thermal mixing. It has been previously employed in studies of flame

acceleration in small-scale experiments.'14 However, when it was used in conjunction with the k-f

model in attempted simulations of flame acceleration experiments, the results were unfavorable.

Strong spatial oscillations of the flow variables and turbulence parameters appeared in the

solution. Although the oscillations were of short wavelength, i.e., alternating on adjacent grid

points, no temporally increasing numerical instability appeared. It appeared that the difficulty was

introduced by trying to resolve phenomena on a length scale incompatible with the grid spacing.

As noted above, this is anticipated when features in the flow configuration are large compared

to the flame thickness. In such cases, in order to achieve a numerically viable description of

the fluid, it is necessary to find some way to define the flame front on a numerical grid which

extends over a large spatial domain. Although there has been some success in the application

of adaptive gridding techniques to this problem, 15- 16 we know of no attempts to apply them to

calculations where a turbulence model must be used to, define combustion rates. The reason for

this is that the resulting turbulent flame profile may not be particularly well-justified theoretically,

so there is no point in devoting a great deal of computational effort to its definition. It is also true

that the geometric complexity of the present experiments would make application of adaptive

gridding very difficult.

For these reasons, We have'developed a combustion model which combines some of the

features of the Magnussen-Hjertager and eddy breakup, models with some artificial definitions

of turbulence parameters in the flame zone for flame propagation at low turbulence levels. This

allows us to resolve the flame thickness over a number of grid points sufficient for good numerical

behavior, and to specify the burn velocity as a realistic function of turbulence level, pressure, and

temperature of the unburned gas. It is important to note that, while we thereby do not compute

the burn velocity from first principles, we are able to achieve fairly good agreement with many

aspects of the experiment with only minimal adjustment of a small number of parameters.

In the next section, we give a brief description of the Conchas-Spray computer code'7

which we have used for the calculations. (The simulations described here have been performed

in a two-dimensional domain, but the model itself is not restricted to two dimensions.) We also

discuss some features of the version of the k-c turbulence model that we have implemented.

into the computer code. In Section 111, we derive the new features of our combustion model and

discuss its capabilities when integrated into the code.

Results of the simulation of an actual experiment in the FLAME facility are presented

in Section IV. The reason for the co ncentration of effort on only one experiment is that only

three experiments have been performed within the particular configuration to which our model

currently applies, and they do not cover a large range of hydrogen concentrations. Both the

experiments and the computer calculations are expensive; this unfortunately limits the data base

and the number of simulations permissible.

A particularly noteworthy result is that there is a period of time during which the computed

flame velocities and pressures agree quite well with the experimental values. 'During that t~ime

perid, flw te vlociiesre pprximately sonic, and it. appers likely that these velocities

are determined primarily by a choking mechanism. .The pressure signals are character ized .by

oscillations determined by the burnout of the individual chambers formed by the obstacles.

Some conclusions and recommendations are offered in the final section of the paper.

4
111. COMPUTER CODE AND TURBULENCE MODEL

Flow Variables and Chemistry

The Conchas-Spray code (modified to include the k-,E turbulence model) solves the following system of equations for the flow variables:

apu) +V(pju) = - V~p±Vlo (2 at

da~pI) + (iu V-v+ u + (P.u [k VT + D)hV~p/~ or + (2) att

where pi is species density, u is velocity, I is specific internal energy, j~and Q~are the rates

of change of pi and energy per unit volume due to chemical reactions, k* is the turbulent kinetic

energy, and E is the rate of dissipation of k. . The index i denotes the individual species. In this

paper four species are considered: H2 , 02, N2, and H120. In the interest of conciseness, the

identification and method of evaluation of many of the rest of the terms in these equations will

only be briefly outlined. For further details, see References 17-19. The total mass density p is

obtained by summing the pi over all species. The temperature T is computed from the energy

I by assuming that the species enthalpies hi are functions only of. T. The pressure p is then

obtained from the ideal gas law. Diffusion of mass and heat are accounted for by the terms

involving the mean diffusion coefficent D and the thermal conductivity kT. The stress tensor

a is given by

p [uVU + VU -- 2(V -u'6 (4) 3 3

where b is the unit tensor, the superscript T denotes the transpose operator, and the viscosity

is

A = AL + AT (5)

where AL is laminar viscosity and MT is turbulent eddy viscosity. The eddy viscosity is obtained

from the k-E turbulence model, and will be discussed later. The laminar viscosity is obtained from

Sutherland's formula with appropriate coefficients, but is actually negligible in our calculations.

The term c in (3) is the rate of dissipation of turbulent kinetic energy per unit mass.- (The

dissipation term describes the conversion of turbulent kinetic energy into heat; hence it is an

5
additive term in the energy equation (3).) Equations (1)-(3) are in Favre-averaged, for,"",

assuming that the turbulence model can appropriately account for the Reynolds stresses via

Equations (4) and (5), and that c can be correctly computed. We are neglecting differential

species diffusion by using an average value of diffusion coefficient D. This is not well-justified

for a hydrogen-air mixture, but in view of the fact that we are using an artificial model to define

the flame zone, this discrepancy does not detract from the results.

The thermal conductivity kT and diffusion coefficient D are related to viscosity it by

kT - CA(6)

D t (7) pSe

where Cp is the sp ecific heat at constant pressure and Pr and Sc are turbulent Prandtl and

Schmidt numbers. We have used Pr = Sc = 0.52. However, because of the aforementioned

flame front modeling, our results are not sensitive to this choice.

The following quantities appearing in the equations remain to be determined: The chemical

contribution to the species conservation equation 'q, and the chemical heat releaseQ*Ths

terms describe processes that are too complex to be evaluated in detail numerically for the

present experimental configuration. Our model -for them will be derived in Section Ill.

In this work, the chemistry is assumed to consist of the global one-step reaction

2H2 + 02 *2H120 (8)

In this case, given the pi's and hi's, all the chemical rates k~and Q r opeeydfndb

specifying only one of the species production rates. (We choose to use *cfor this purpose.)

At this point, it is convenient to define one more variable for later use. Models which

assume that combustion is limited by the rate of thermal mixing usually make use of a progress

variable, or something related to such a quantity. Our progress variable is defined as

where Y,, is the mass fraction of the product species and YPOis the mass fraction of product

that would result if the locally extant mixture was allowed to react to completion. (Subscripts f

and p will be used for fuel and product henceforth. In our case these correspond to H2 and

H120.) Then

Pp + VPPP 700 = If Mff

p

where the vi are stoichiometric coefficients (both equal to 2 in our one-step H2-air reaction-see

Equation (8)), and the mi are the corresponding molecular weights. Hence,

6
CPp _ p (9) C=pp + pm~pp flp + (9)ln

where flf and flp are the number densities of the fuel and product species. The quantity C

always satisfies 0 < C < 1. (Note that the denominators in Equation (9) are the densities

of product species instantaneously available under conditions of complete reaction of a lean

mixture (at constant volume), and C = 1 when the fuel is completely consumed. This equation

can be .used for a fuel-rich mixture; if the subscript f is replaced by o (denoting the oxidizer),

C = 1 when the oxidizer is consumed.)

The computer code used used for our calculations is an adaptation of the Conchas-Spray

code developed by Los Alamos National Laboratory. It provides a solution to Equations (1)-

(3) in two dimensions; the present calculations are performed in rectangular coordinates. The

numerical algorithm employed in the code is a finite-difference approximation derived from an

integral formulation of the conservation laws. For details see References 17 and 20.

All the computations employ an Eulerian method on a grid with square cells (see Figure 3). The thickness of the obstacles is taken to be two grid spacings for numerical simplicity.

(The experimental obstacles were constructed of 13 mm (1/2 inch) plywood, and were therefore

considerably thinner. However, the difference has negligible effect on our results.) The ,;"axis

lies in the midplane of the facility, which is assumed to be a symmetry plane. Hence, the grid on

which the computations were carried out encompassed only half the area of the facility. Except where noted otherwise, the grid consists of 13 points in the (transverse) x-direction by 397 grid

points in the (axial) z-direction. This encompasses sixteen and one-half chambers formed in the

channel by the obstacles used in the experiments.. This is a very coarse grid; it is not adequate

to resolve recirculating flows behind the obstacles with precision. Furthermore, errors will be

incurred in the computation of the shock waves which form when the gas flows past obstacles

at supersonic speeds. We have investigated the numerical accuracy to the extent possible by refining the grid spacing by a factor of 2/3. The results (discussed in Section IV) show some

differences; nonetheless, the coarse grid suffices for calculations which provide insight into the

processes involved.

The k-E Turbulence Model

The turbulence model which we have used is the k-c model as described in References

9 and 10. Although Reference 9 specifical~ly addresses the question of compressible flows, it

should be noted that turbulence models for compressible flows are not well-developed. Hence,

the, choice of the standard k-f model in this situation cannot be regarded as definitive.

The model consists of the following transport equations for the evaluation of turbulent

kinetic energy k and dissipation rate f.

at

(P)+ V - (pCU) = C1, (Cr:V U) + V.- (

1-Vf) - 2P 1) atpE kc or, .

7
,396 zones

a a

X24 zones

!a 12 zones

Figure 3. Problem domain and finite-difference grid used in most of the calculations. The full

16! chambers into which the facility is divided is shown on the left. The grid is too fine to be 2

resolved on the scale used there. An expanded view of 3 chambers is given on the right.

8
The terms on the right sides of these equations represent (from left to right) the production,

diffusion, and dissipation of the quantities k and c. Strictly speaking, Favre averaging requires

the inclusion of a term in (11) involving the correlation between velocity fluctuations and the

pressure.9"13 Following a practice often used'3 , we omit this term to avoid undue complexity in

a turbulence model which is still under development for compressible flows. The coefficients

C1, and C2 , are modeling constants, and or, is a Schmidt number. These. constants have been

tuned'0 for agreement with certain experiments, with optimum values C1, = 1.92, C2, = 1.44,

and or, = 1.3. (The experiments used for tuning are much simpler than those addressed here.

Furthermore, they involve incompressible flows. Hence, the use of these. constants should not

be. regarded as definitive.)

The turbulent part of the viscosity (the eddy viscosity) is given by

AtT CAP (12)

The coefficient C.~ is set equal to the usual value of 0.09 (see, e.g., Reference 10).

Since k and c are defined to be nonnegative, it is important that the numerical scheme.

preserve this property. The production and dissipation terms in the k and c equations require

special consideration in this regard. Since turbulent regimes are sometimes characterized by

the fact that production approximately equals dissipation, care must be taken to ensure that

numerical errors in either term do not result in negative k or E. A special numerical scheme was

devised to. alleviate such difficulties. It is given in Appendix A.

In the Conchas-Spray code, the boundary conditions on tangential velocity are obtained by

implementing the law of the wall for turbulent boundary layers.'17

' 21 ,

22 To this end, the tangential

velocity at walls is not zero, but is allowed to vary to account for fluid momentum in the wall grid

cells, which are. in the boundary layer. We have essentially retained the numerical boundary

condition as in the original code except that we restrict the velocity in a cell at a wall to be less

in magnitude than the tangential velocity in the adjacent cell one grid point in from the wall. This

prevents occasional nonphysical profiles from occurring.

The, boundary conditions on k and c are likewise obtained from consideration of the lawof-the-wall shear layer. The wall shear stress r is obtained by assuming that the tangential

velocities at grid points one cell away from the wall are in the law-of-the-wall boundary layer.

This allows one to solve a simple equation for the shear stress.'17 Assuming that the Reynolds

stresses are nearly constant, it can be shown10 that the turbulent kinetic energy in the wall cell

may be approximated by

VC'

where u,~ is the shear velocity

u*= \r1p

9
This provides the boundary condition on k. The associated boundary condition on f is

3

where y is the distance to the wall and Kc 0.4 is the von Kirm~n constant. This Il/y behavior

of E can give rise to some numerical inaccuracy, especially in the diffusive flux of C. This can

be. alleviated by implementing a special differentiation formula near the wall. Details are given

in Appendix B.

11ll1. THE COMBUSTION MODEL

Properties of Turbulent Flames and Numerical Flames

As noted in the Introduction, simulation of combustion processes in large experiments

places extreme demands on the capabilities of even the largest and fastest computers. In

the present case, it is necessary to give up the idea of resolving the flame thickness on the

computational grid. A possible exception arises in situations in which the turbulent length scale

is very large and the turbulent flame thickness, defined as the thickness of a region where

burned and unburned mixture both exist, extends over several grid points. But for satisfactory

numerical behavior we must usually thicken the flame artificially. To illustrate this point, a plot

of the laminar, turbulent, and computational flame profiles are given in Figure 4. Turbulence

results in rapidly fluctuating wrinkles and spirals of interpenetrating burned and unburned gas,

as shown in the upper sketch in the figure. The turbulent flame thickness 6b. is defined by the

region where the composition is intermittent due to the fluctuations, and the average value of C

is greater than zero, but less than unity. For good numerical behavior, the computational flame

thickness b, must be equal to a few grid spacings. (The size of one grid interval is indicated by

A~x. The scale in the figure is not necessarily accurate; the various lengths actually satisfy

be ': AsX > 6b' > 6L 4

Artificial flame thickening has been done previously by a method known as the /

transformation. 23' 24 This involves changing the length scales in the problem by artificially increasing the thermal diffusivity and diffusion coefficient and decreasing the combustion rate.

We have retained this essential idea. However, the original #3 - transformation established the

scaling through the use of the local temperature gradient in the problem. Numerical experimentation indicated that this did not work well with mixing-limited combustion models such as

the Magnussen-Hjertager model. Hence, it was decided to achieve the flame thickening in a

different way.

Given that the flame will be somewhat artificial, it is necessary to first determine what

properties it should have. In particular, it is desirable that the burn velocity be realistic. The

10
Burned BurnedUnburned

-- Lamninar --- Turbulent - - - Computational

I

H-- AX P.I

Figure 4. Schematic illustration of laminar flame thickness, turbulent flame thickness, and computational flame thickness. Above: A slice through a turbulent flame, plotted in (x,y) coordinates.

The thickness of the dark lines is the laminar flame thickness'bL. Below: A plot of the progress

Variable C (defined in Equation (9)), for y equal to the constant value defined by the dashed

line in the upper figure.

11
following properties of turbulent flames have been extracted from a list given in Reference 13

(SL and ST refer to laminar and turbulent burn velocities, respectively):

1. ST > SL (almost always)

2. ST increases with u' Vk. Often, the increase is linear, and

ST ASL +Bu'

provides as good a correlation as anything else*.

3. Sometimes ST is insensitive to the scale of approach flow turbulence.

4. In open flames, the variation of ST with composition is much the same as for SL.

5. In flames in ducts, ST is larger than for open flames and is sometimes insensitive to

approach flow turbulence and composition.

Our model represents an attempt to adhere to these pr operties, except for the latter clause

in. Number 5. It does not necessarily apply to our experimental configuration, as it was inferred

from experiments on steady-state flames stabilized in small ducts.

To see how to embed the above properties 1 through 4 in a combustion model, we note

that, for the case of mixing-limited combustion models such as Magnussen-Hjertager or eddy

breakup, simple flat flames will be expected to have the following properties (see, for example,

the discussion of laminar deflagrations in Section 5.1.2 of Reference 19):

The turbulent burn velocity will scale as

S X/tT (1.3)

where X is the thermal diffusivity (assumed proportional to the diffusion coefficient-see Equations

(6) and (7)) and tT is the turbulent time scale. (This is commensurate with the assumption of

burning limited by thermal mixing, where the rate of combustion is proportional to 1 ItT).

The turbulent flame thickness will scale as

b~ XtT (14)

From the k-c turbulence model,

x- X C~V4P

and

*Note that u' is *defined here as the average turbulent velocity, and not as the instantaneous

velocityfluctuation.

12
L (5 tT (5

where

L k3/2/C (1.6)

is the integral length scale of turbulence. From (13) and (14), then

ST oc Vk u (17).

b cxL (8

In other words, if the turbule nce intensity v~k and the length scale L are given for the gas, one

.can achieve a model which complies with the linear part, of Property 2 and with Property 3.

Development of a Model With Appropriate Properties

In applying a mixing-limited combustion model to a computation such as this, it must be

recognized that difficulties may occur when starting up a problem when the turbulence level is

low. Initially, one expects the flame to be laminar, but in our experimental configuration this

will be quickly superseded by a flame which exhibits some form of turbulence, e.g., a'wrinkled

laminar flame.' 3"19 Data from a variety of experiments in the FLAME facility suggest that the

hydrogen-air flames studied in the facility achieve a relatively constant burn velocity Which is

greater than laminar and which persists until large-scale instabilities and/or obstacles cause the

flamhe to accelerate.2

' 3

',25

The approach that we have developed is to endow the gas with minimum values of k

and L, denoted k.. and L,, which are used for the computation of thermal conductivity and

diffusivity only. In that way, according to Equations (17) and (18), there will be minimum values

of burn velocity and flame thickness. The minimum burn velocity provides a term like ASL in

property 2, although not in such a way that the equation there is exactly satisfied. The scheme

which we employ is to use km, and L.. in the formulas for kT and D unless the true turbulence

intensi ,ty k or length scale L predicted by the k-c model exceeds those threshold values. Then

we use k and/or L. It is crucial to note that only the true values k and L (or, equivalently, E) are

used in the viscosity. In particular, they are used in the transport equations (10) and (11) for k

and E themselves, and in (12) for ALT, so that km, and Lm, affect the generation of turbulence only

indirectly., (it may appear that this indirect effect is a very strong one, since the gas velocities

appear Iin the turbulence. production terms in (10) and (11). However, if km and Lm' are chosen

correctly, these gas velocities would be essentially the same as those. occurring in a perfectly resolved flame.)

Now, how should km, and L. be determined? It would be desirable that km correspond

to something physically realistic, because the point at which k exceeds km should represent

some real point in parameter space at w hich the flame starts to propagate more rapidly.

13
We have extracted a reasonable choice from the work of Libby, Bray, and Moss (LBM) 26

In the limit of large turbulent Reynolds number RT pu'L/II, they obtain a formula.

ST SL + 1.14 Vik (19)

This formula is derived by making certain modeling assumptions regarding physical processes

in flames. It is not intended to be universally valid, even in the limit RT --* 00. In fact, References

13 and 26 demonstrate clearly that no universal formula exists. But we choose to use this as a

starting point. Given that we can extract an initial turbulent burn velocity STj from experimental

data, we obtain a provisional estimate for km by inverting (19):

km (T S) (20)

We emphasize that this is only an estimate. As discussed- below, it is impossible to satisfy all

the desired flame properties and still adhere to this value of k.m.

We now give our formula for the rate of consumption of fuel. (As noted in the discussion

following Equation (8), for a one-step reaction, specification of the. rate for only one species

suffices to specify all chemical rates..) Except for a dependence on pressure and unburned gas

temperature to be discussed. later, the combustion rate is

=f -Acpf (T)C( CO/T (21)

where A, is a dimensionless constant, and

T,<T

f(T)={T-TO 0-< <T

0 ~T <T

where T, is adiabatic flame temperature and T, is the initial gas temperature. Note that the

factor C(1 - C) is similar to that appearing in approximations to the eddy breakup model (see

Equation-(4.72) of Reference 13) and to the minimum normalized mass fraction appearing in

the. Magnussen-Hjertager model'". The function f (T) serves a spec *ial numerical purpose. It

does not really represent any physical process, except that it provides a qualitative. increase

in combustion rate with temperature up to the adiabatic flame temperature. It is inserted in

Equation (2 1) so that the flame profile will be determined by thermal diffusion. in such a way as

to obey Equations (13) and (14), with the result that (17) and (18) are satisfied. Lacking the factor

f (T), Equation (21) results in flame zones that tend to diffuse too far ahead of and behind the

flame, and which do not exhibit the functional behavior specified in (13) and (14).. This spreading

-is predictable on the basis of comparison with solutions to the Kolmogoroff-Petrovsky-Piscounoff

(KPIP) equation21

7 29.

14
At this point we have three adjustable parameters: kin,, A, and Lm. They are fixed by

simultaneously specifyingf the following three flame properties: The initial turbulent burn velocity

5T:, the derivative SIT = dST/dVk- for adiabatic flames, and an appropriate numerical flame

thickness. Our initial choice for dST/dv'ik is 1.14, in accordance with (19). A variation on this

value will be discussed in Section IV. (The actual values of kin, A,, and Lm are determined by

a djusting them in a series of one-dimensional flat flame calculations.) In most of our calculations,

we have attempted to maintain a 10-90% thickness of 4 grid points., As noted,.above, our model

does not precisely satisfy Equation (19). It does, however, provide values for the initial burn

velocity and slope of the ST vs kkcurve which are physically well-justified. (See Figure 5.)

Dependence on Pressure and Pret lame Tlemperatuire

Finally, we describe the-way in which the combustion model is augmented to include a

dependence on pressure and the temperature of the unburned gas through which the flame

is propagating. As the flame travels down the obstructed channel of the FLAME facility, the

unburned gas in front of it is compressed due to the overall drag forces which build up. This

compression results in an increase in pressure. and temperature. Since the laminar, burn velocity

of a hydrogen-air mixture depends on the pressure and temperature of the unburned gas, the

turbulent burn velocity will also exhibit such a dependence. In the spirit of property 4 which is

desired for the combustion model, we simply choose to impose the same functional dependence

Of ST' on p and unburned gas temperature T, as is expected for the laminar burn velocity SL.

WarnatZ30 has shown that in certain regimes SL varies according to

SL C P

t T.~

Consideration -of Equations (13) and (14) indicates that, all else remaining constant, ST Will exhibit the same behavior if the burn rate, thermal diff usivity, and diff usion coeff icient are multiplied

by a factor

h (p, T,) =(22)

where p, is the initial pressure and T, is the unburned gas temperature. We have used the

values m = 0.2,n = 1.64 as being representative of Warnatz' results; this choice is by no

means clear-cut, however.30

We have employed Tug the unburned gas temperature, to alter the burn rate and diffu'sivities. But ', X, and D vary thro .ugh the flame. How do we know what 72,. is at points inside the

flam? Ths isdetermined by extrapolating back to an unburned mixture at the same pressure. (It is assumed that the pressure does not vary significantly through the flame.) We write

jcp(T)dT = A (23)

T. MPy

15
- Model (S'T = 1.14)

------- Variation in

model (S'T

- - -LBM theory

= 0.82).

ST

STSL9

kmV

VVk

Figure Turbulnt 5 burnveoiyv uruecinestfoth ib-a-M smd12(hiFigur 5. Tubietbrnvlciyv turbulence intensity ro the forls Libby-Brnueriay-os moelechanimaalbein-th

present case, and not from any suggestion that the physical situation is better represented.

The dashed line represents a variation in the model in which the combustion rate is reduced

(see discussion in Section IV). The scales in the figure represent a realistic range of values

encountered in the calculations; typical values of N~ in the choking regime are about 30 times

16,
where cp(T) is the specific heat of the mixture and AýH is the molar heat of reaction (57..8,

Kcal/mole for H2 0). We have obtained a linear fit to the specific heat from data available in the

CHEMVKIN 31 code. This approximation has the form

CP = PO +c'(T - T,) (24)

When (24) is substituted into Equation'(23), there results a quadratic equation for T.; the solution

is substituted into Equation (22) to obtain the function h(p,T,).

Summary ot the Model

This completes the description of the combustion model. The implementation of the

material presented in the previous section requires a redefinition of the thermal conductivity and

diffusion coefficient and an alteration of the formul a for combustion rate given in equation (21)

For convenience, the appr opr iate formulas are collected in the following:

kT =CpMIt('

Pr

Pf =-Acp f (T)h(p,T,)G(1 C)/t'T (21')

where

A1 A L + (5'

A'T C=p Vk L, (12')

and k, and L, are effective values of k and L, defined as

ke =maxf k,km} (25)

Le =max{L,Lm} (26)

The model has been tested by running one-dimensional calculations of adiabatic flames

.at constant pressure. For this purpose, the k-c model was turned off so that km and Lm would

17
determine the flame properties. In-such calculations the flame achieves a self-similar profile

propagating at an approximately constant velocity. As indicated in the foregoing discussions, the

desired behavior of the flame velocity and the computational flame thickness b, vary according

to

ST oc \/kp' T. (27)

b, CcL (28)

Furthermore, ST and b, should be independent of the grid size. Table I gives some results obtained by varying the parameters around a representative baseline set. (Anticipating

the simulations to be discussed in Section IV, the baseline calculation corresponds to a 14.5%

hydrogen-air burn at Atmospheric pressure and 3000 K with k,,in, Lm, and A, chosen appropriately for one of the calculations.) The agreement between the computational values and the

desired values is quite good, except possibly for the pressure and temperature dependences.

The reason for the discrepancy in the case of the temperature variation is that there is some

inaccuracy in the determination of T,, via Equations (23) and (24). Since there is some arbitrariness in the choices for m and n which are used in Equation (22), these discrepancies will be

accepted for the remainder of this paper. For the most part, the model achieves its goals for

the simple flat (one-dimensional) flame. In particular, it may be noted that the first three lines in

the table verify that the slope dST/dvlk- is approximately 1.14, as required by Equation (19).

In summary, the characteristics of the Model are as, follows: It is designed to provide for

the consumption of fuel at a rate consistent with certain assumptions as to the phenomenological behavior of turbulent flame propagation. It permits the computation of the behavior of~a

compressible fluid ahead of and behind the flame front; it thereby accounts for the interaction

of turbulence with the flame insofar as the assumptions relating the turbulent time scale to the

combustion rate are correct. It does not permit a computation of burn velocity from a knowledge

of only the state of the gas and the turbulent flow structure; it requires prior assumptions as to

what the burn velocity would be in a given situation.

It should also be noted that, in the form presented here, the fluid model will exhibit a nonphysical heat conduction and particle diffusion everywhere because of the universal application

of the (artificially high) threshold values of k and L in Equation (12'). For the present application,

errors due to this effect are negligible. To see this, consider the energy conservation equation

(3). We wish to show that the artificial heat conduction is negligible on the time scales typical of

changes in the specific energy. Suppose that the typical time scale tt is the time for the initial

flame to propagate the width of the channel:

t tVfO

where w'is channel width and v1 j, is the initial flame speed. If w is used as the. thermal diffusion

length scale, the ratio of the diffusion time t~ to tt is

18
Table 1.

Results of tests of the combustion model. The quantity Aýx is the

grid spacing. The burn velocity and flame thickness resulting from the

model are ST and b4; the values of these quantities expected on the

basis of simple theory are denoted S~ and b,". The constant A, (see

Equation (21) is equal to 30.8. The exponents m and n (see Equation

(22)) are equal to 0.2 and 1.64.

k.(J/lkg) L.,(cm) p(atm) T.(OK) Ax(cm) S.e(mls) ST(m/s). 6c(cm) 6. (cm) Comments

4.7 68 1 300 7.62 - 2.5 -16 Baseline Case

19.0 68 1 300 7.62 5.0 4.9 16 16 km 4k~ae

1.2 68 1 300 7.62 1.25 1.24 16 16 km I 'bs

4.17 136. 1 300 7.62 2.5 2.4 32 34 Lmn2L!:86

4.7 68 7.6 300 7.62 3.7 4.4 16 25 (p/pbas)m =1.5

4.7 68 1 458 7.62 5.0 4.4 16 19 (T/Tbase)n -2

4.7 68 1 300 3.81 2.5 2.4 16 16 LX- Azbdae 2

4.7 68 1 300 15.24 2.5 2.6 16 18 1Ax 2. Axa.e

19
t W 2

tt Xtt

vf~wPr

CovfkmLm

For typical values vf, 16 m/s, w =1.83 m, Pr 0.52, Cm 0.09j,t km,= 4.7 J/kg,

Lm =0.68 m, this ratio is 115. Hence, the diffusion time is long compared to the time during

which other processes are occurring, and the artificial diffusion is negligible (except in the flame,

.of course). As time elapses and the flame accelerates, tt decreases and the relative error

incurred becomes smaller. Inspection of the dispersion relation for acoustic waves shows t hat

the attenuation of such waves. is similarly negligible over a propagation -distance of a channel

width. The attenuation of acoustic waves that propagate the entire length of the FLAME Facility

might not be negligible, but they would not be, expected to contribute significantly to the results

of the present study.

The effects of this nonphysical diffusion could probably be eliminated entirely if necessary.

All that would be required is to remove the thresholds on k and L away from the fla me. This,

would require some care, as it would be necessary to ensure that the resulting decrease in

diffusivities -was kept far enough away from the flame -zone that the flame propagation velocity

was, not affected.

Finally, we note that no attempt has been made in this work to account for the effects of

flame stretch and flame curvature in detail. (A general discussion is given in, e.g., Reference

19.) It is impossible to compute such effects in this case, since the flame is artificially thickened,

but one might attempt to model them (see, e.g., the model for flame quenching in Reference 4

32).

IV. COMPARISON OF COMPUTATIONAL RESULTS WITH EXPERIMENT

Experimental Considerations

A total of 9 experiments in which obstacles were present have been performed in the

FLAME facility. Table 11 gives a brief description of these tests and their results.33'

34 in Tests F21 through F-23, the facility was unvented; in F-24 through F-29, one-half of the top plates were

removed, resulting in 50% venting. In the present report, we are considering only the unvented

experiments. In experiment F-21, failure of the mixing fans resulted in uncertain hydrogen

concentrations. Hence, only experiments F-22 and F-23 are available to us for comparison with

the type of computations with which we are concerned.

The only difference in the experimental arrangement in F-22 and F-23 is that the initial

hydrogen concentrations were:15% and 14.5%, respectively. In experiment F-23, the mixture

appeared to detonate;. in F-24 it did not. Since turbulent flows and conditions leading to detonation are random in nature, the 0.5% difference in concentrations in these two experiments

20
Table 11.

Summary of FLAME tests with obstacles. (From Reference 34.) A total of 16 pairs of

obstacles with a blockage ratio of. 1/3 were used in the experiments. (See Figures 1 and 2.)

Peak Equivalent

Test Top Hydrogen Peak Planar

Number Venting Mole Overpressure, flame Speed

Area,% Fraction,% atm rn/s

F-21 0 10-15a 6.5 580

F-22 0 15.0 31.0 700

F-23 0 14.5 12.0 540

F-24 50 15.5 b 46

F-25 50 19.7 16.5 *890

F-26 50 285 19.7 16

F-27 50 13.1 .09? 15?

F-28 50 14.9 .09? 33.4

F-29 50 18.5 .22? 130

.a - mixing fans did not operate.

b - signal in noise level.

-analysis not complete.

21
is not necessarily significant. In other words, it is not necessarily true that a 15% concentration

would always lead to a deflagration-to-detonation transition (DDT), while a 14.5% concentration

would not. Our computations are not capable of the spatial and temporal resolution or the

detailed physics and chemistry required to provide a true simulation of DDT. Nonetheless;,it will

be seen that we observe a "numerical" transition. Hence, we will concentrate on a comparison

with experiment F-23, with a brief observation on detonation with a view to experiment F-22.

In the FLAME facility, flame arrival times are measured by detecting temperature increases

at 5 thermocouple "rakes" positioned vertically on the channel midplane. (See, Figure 2.) In the

present set of experiments, each, rake has thermocouples at 7 vertical locations. From this

arrival time data, it is possible to construct rough vertical profiles of the flame as it propagates

down the channel .

33 The arrival time data for experiments F-22 and F-23 and the corresponding

flame front profiles are given in Figures 6-9. The flame velocity obtained from the slope of

the line connecting the last two da ta points in Figure 6 is only 700 m/s. This is below the

Chapman-Jouguet (CJ) detonation velocity of 1520 m/s, but it is believed that the discrepancy

is due to the poor resolution available from the thermocouple data. Presence of a detonation

was inferred from visual observation. Note also the rapid propagation of the'flat profiles which

begins sometime between 300 ms and 320 ms in Figure 7. This is indicative of a detonation.

In addition to the thermocouples, there are pressure transducers located at various sites

on the side walls. Dat a from these transducers will be compared with the computational results

below.

Comparison of Flame Trajectories and Pressure Histories

The first calculation carried out to simulate Experiment F-23 was made with the coefficients

in the combustion model selected precisely as described in Section 111. (it will be referred to as

the baseline calculation.) Parameters pertinent to the test configuration and the computation

are collected in Table Ill. Ignition in these computations is accomplished by forcing a controlled

burn over a small region at the closed end of the channel. This ignition process ends at a time

of 50 ins. Some care must be taken in implementing ignition in these simulations; a discussion

of this point is given later in this section. The initial computational burn veloc~ity ST1 was chosen

as 2.5 m/s so that the initial computational flame speed would correspond to the experimental

value of 20 m/s. It is to be emphasized that, for this calculation, this is the only parameter

adjusted to agree with experiment. Note that the flame speed expected for an adiabatic flat

flame with a constant burn velocity in an unobstructed tube would be

VT = rTSi

where r is the volume expansion ratio. For a 14.5% hydrogen-air mixture, r is 4.45. Hence,

according to the above formula, VT would be 11 m/s. The reason that the experimental flame

speed is greater than this is that the flame is not flat, but is stretched out, primarily by the

obstacles. As an additional point of reference, it may be noted (see Table Ill) that the laminar

burn velocity at this concentration is 0.9 m/s. The existence of a wrinkled laminar flame or a

weakly turbulent flame with burn velocity equal to, say, 2.5 times the laminar burn velocity, is

consistent with other experimental results and theoretical considerations13 9

22
350-

300-

TC Ee.

250 El 0.076.

IL 0.381

LL 0.762 2i 20 .1 1.143

0 1.524

EB 1.905

2.286

150-

0 5 10 15 20 25 30

DISTANCE FROM IGNITION END, m

Figure 6. Flame arrival time data for Experiment F-22. (This and-the following 3 figures were

provided by Reference 3ý4.) The elevations of the thermocouples are given in 'the key. The

distances of the thermocouple rakes from the closed end of the channel are: 1.0 m, 3.8 m, 11.6

m, 19.8 m, and 37.8 mn (see the locations of the sets of data points).

23
TO SCALE

2

E

2f 0

1=

wj -J

1.5

0.5

0 510 15 20 25 30

DISTANCE FROM IGNITION. END, m

Figure 7. Flame front profiles for Experiment F-22. A side view of the facility is shown (the section.

is through the midplane, where the thermocouple rakes are located). The upper figure is to scale.

The lower figure is expanded in the vertical direction, showing the times in milliseconds to which

the contours correspond. The data points indicate the vertical positions of the thermocouples

(but not their axial locations-these profiles are obtained by interpolation).

A

24
NUREG/CR-4855

SAN D87-8203

R3

DEVELOPMENT AND APPLICATION OF

A COMPUTER MODEL FOR LARGE-SCALE

FLAME ACCELERATION EXPERIMENTS

K. D. Marx

July 1987

Sandia National Laboratories

Livermore, CA 94550

Operated by - Sandia Corporation

for the

U. S. Department of Energy

Prepared for

Division of Reactor Accident Analysis

Office of Nuclear Regulatory Research

U. S. Nuclear Regulatory Commission

Washington, DC 20555

Under Memorandum of Understanding DOE 40-550-75

NRC FIN No. A1246
A
400

350

co)

E

U0

w

.300

250

200

150

100

0 5 10 15 20 25 30

DISTANCE FROM IGNITION END, in.

Figure 8. Flame arrival times for Experiment F-23. (See Figure Caption, 6 for explanation.)

25
TO SCALE

2

E

2f1.5 0

LU 1

-j

0.5

0 5-' 10 15 20 25 30

DISTANCE FROM IGNITION END, m

Figure 9. Flame front profiles for Experiment F-23. (See Figure Caption 7 for explanation.) No

detonation occurred in this case, although in general the flame propagation rates appear to be

only slightly slower then for F-22.

26
Table Ill.

Summary of parameters pertinent to Experiment F-23 and the numerical simulation.

Hydrogen mole fraction

Initial density

Initial temperature T,

Initial pressure

Initial sound speed

Adiabatic flame temperature

Volume expansion ratio r for adiabatic flame

Laminar burn velocity SL

Isobaric sound speed

Isochoric temperature

Isochoric pressure

Iscohoric sound speed

Chapman-Jouguet detonation velocity

Chapman-Jouguet pressure,

Initial experimental flame speed

Combustion rate coefficient A,

Threshold turbulent kinetic energy km,

Threshold turbulent length scale L..

Pressure exponent m

Temperature exponent n

.Prandtl number Pr

Schmidt. number Sc

Initial turbulent burn velocity STi

Initial computational flame speed

0.145

1.0 14 kg/rn3

3000K

1 atm

374 m/s

1438 OK

4.45.

0.9 m/s

758 m/s

1757 OK

5.43 atm

834 m/s

1500 m/s

11 atm

20 m/s

30.8

4.74 J/kg

0.68 m

0.2

1.64

0.52

0.52

2.5 m/s

20 m/s

27
It is interesting to note that the burn velocity of 2.5 rn/s for a hydrogen concentration of

14.5% is between 1/2 to 2/3 of the value computed independently from an analysis of experiments without obstacles"5. Since non-planarity of the flame was ignored in obtaining the latter

results, this lower value is reasonably consistent.

IA comparison of the computed and experimental flame trajectories of one calculation with

the corresponding experimental data taken from Experiment F-23 (14.5% hydrogen) is given

in Figure 10. Considering the size of the experiment and the complexity of the computational

problem, the agreement is quite good. To some degree, the discrepancies in the computed

flame trajectory must be viewed in light of the fact that the flame velocity is extremely sensitive

to the burn velocity. (Because of the geometrical factors and the feedback loop involving gas

flow velocities, turbulence intensities, and combustion rates, the dependence is nonlinear.) This

sensitivity. notwithstanding, the computation underestimates the rapid acceleration of the burn

rate in the region from 200 to 275 milliseconds. The reasons for this are not known. Precise

comparisons with experiment are not possible, as no measurements of turbulence levels are

available and it is not possible to instrument the experiments sufficiently to make detailed spatial

plots of such flow parameters as flame position, pressure, and temperature.

At a time of about 300 ms into the calculation, however, the computation is accelerating

rapidly. Eventually, the computational flame trajectory essentially parallels that of the experiment

while the flame traverses through axial positions from about 17 m through 22 m.

Figures 11-14 give a comparison of the computed and experimental pressures at the

locations indicated in Figure 10. Data from two pressure transducers at axial positions 18.6

m, 21.5 m, and 25.6 m are shown. Table IV gives complete information on the locations of all

transducers pertinent to this paper. The computed pressures are taken from the finite-difference

zone next to the channel wall. (Note: All of our calculations terminated when their flames

reached the end of the channel. That is why the computed pressure histories end abruptly at

3.20 ms in Figures 11-14.)

The agreement between computation and experiment at 14.6 m is rather poor, as seen

in Figure 11. This simply illustrates the point that the computed combustion rate is significantly

slower than that occurring in the experiment at that location. (Compare the flame speeds, i.e.,

the slopes of the trajectories in Figure 10.) It is important to realize that the agreement with the

experimental pressures in this figure could easily be improved by adjusting the parameters kin,

Lm, and A, until the flame speed in the region around z = 14.6 m approximately coincided with

the experimental flame speed, but then the agreement would not be as good elsewhere (see

below).

A significant improvement is apparent in Figure 12. The computed pressures agree reasonably well with those measured experimentally. The reason for this is that the structure of

the pressure history is characterized by the timing of the burnout of the gas in the successive

chambers. The fact that the flame velocities approximately agree means that the chambers are

burning out at about the correct rate-. It should be noted that the flame velocity of about 540 m/s

is greater than the sound velocity in the unburned gas, but less than that in the burned gas. This

is to be expected if the flow rates are reduced due to choking6. Note also that the experimental

pressure signals from different transducers at the same axial location are somewhat different.

28
40

-wi--Experiment F-23 ......... Calculation

.....-- Location of pressure transducers

30

E

0.

.0 0

1u4

0-

0 100 200 300 400

Time (ins)

Figure 10. Comparison of computed and experimental flame trajectories for Experiment F-23.

The experimental data has been shifted 100 ms to the left (relative to Figures 5-9) to facilitate

comparison. (In view of the lack of a detailed simulation of ignition, the relative origins, of the

computational and experimental time scales is somewhat arbitrary.) The calculation used the

baseline parameters as indicated in the text and the first line of Table 1. The horizontal dotted

lines show the location of the pressure transducers (at approximately 14.6 m, 18.6 m, 21.5 m,

and 25.6 m down the channel) which provided the data discussed in this paper. (The complete

array of experimental pressure data is not restricted to these locations.)

29
8

-Experiment (PKU4)

........ Calculation 6

,On

E)

4

2-

0--

280 290 300 310 320 330

Time (ins)

Figure 11. Pressure histories at 14.6 m from the baseline calculation and from transducer PKU4

in the experiment. (See Figure 10 and Table IV for a precise description of transducer. location.)'

In this and all other pressure-time. histories, the experimental signal has been shifted in time so

that the times at which the flame passes the pressure transducer agrees with that obtained from

the computation.

30
8

6

IL

Experiment (PKU5).

- -- -- Experiment (PKU1 0) .......... Calculation

4

gsa,

I ~I *~*

g ~ 14

a',

I

a

2

0

300 305 ' 310 315 320 325

Time (ins)

Figure 12. Pressure histories at 18.6 mn.

31
12

10

Experiment (PKU6)

-- -- - -- Experiment (PKU9)

.......... calculation

(on wall)

- ' -Calculation

(U (on axis)

0

1..

0

(I)

0

1~ 0.

8

6

4

2

0

9

I

'a

300 305 310 315 320 325

Time (ins)

Figure 13. Pressure, histories at 21.5 m. In this and the following figure, pressures computed

on the axis of the facility are included for comparison (see text).

32
50

40

00t

EU

30

20

z = 25.6 m

Experiment (PKU7)

--------------- Experiment (PKU8)

............... " ..Calculation (on wall)

Calculation (on axis)

315

4

I'

320 .325 330 335

Time (ins)

Figure 14. Pressure histories at 25.6 mn.

33
Table IV.

Location of pressure transducers in Experiment F-23. (From Reference

34.) The coordinates x, y, and z correspond to spanwise, vertical, and axial

directions. (Positive and negative values of x indicate transducers mounted on

opposite sides of the channel.) All distances are in meters.

Transducer xyz

PKU4 -0.914 0.30 14.6

.PKU5 -0.914 0.31 18.6

*PKU6 0.914 2.12 21.5

PKU7 0.914 0.29 25.5

PKUB -0.914 2.12 25.6

PKU-9 -0.914 2.11 21.4

PKUIO 0.914 0.21 18.7

34
This illustrates the complexity of the problem; these signals may differ because of different proximity to the floor or ceiling, or simply because of random fluctuations in the combustion and

Pressure waves.

it may be noted that the computed pressure signals do not show as much high-frequency

structure as the experimental data. This is undoubtedly due to some numerical dissipation exhibited by the computational algorithm. This will be discussed below with regard to a computation

performed on a finer grid.

The order of magnitude of the pressure pulses which occur when the flow is choked,

or nearly so, can be estimated as follows: The peak burn rate in the corners in front of the

obstacles is very high. Hence, it is reasonable to suppose that the burn there is essentially

isochoric. As noted in Table Ill, the isochoric pressure for 14.5% H2 at 3000 K is 5.4 atm (65

psi). This is seen to be roughly verified both experimentally and computationally in Figures 12

and 13. There is naturally some spreading due to. dispersion, especially at locations far from

the flame,,and this reduces the amplitude of the pressure signals.

The agreement between computational and experimental pressures is not as good at the

21.5 m station (see Figure 13). The reason for this is that the numerical flame is burning more

rapidly than was the case in the experiment; in fact, at this point the numerical flame is very

close to a transition to detonation. It is interesting to note that for the case shown in Figure

13, the pressure computed on the axis of the facility agrees more closely with experiment than

that computed in the zone next to the wall, where the experimental pressure is measured. The

reason for this is believed to be that the computed combustion rate at the wall is too high, and

that this is causing the transition to detonation.

Figure 14 shows the pressures at 25.6 mn down the channel. Detonation has occurred in

the computation. As noted above, no detonation was observed in the experiment, but the peak

pressure of 13 atm (12 atm overpressure, as indicated in Table 11) exceeds the Chapman-Jouguet

pressure of 11 atm. The peak computed pressure of 48 atm far exceeds the CJ pressure as

well. This may be due in part to numerical overshoot, but it is also true that in the computation,

symmetric waves are reflecting off the side walls. Since the pressure in a reflected wave adds

in phase and multiple reflected waves can occur, there is a strong mechanism leading to the

computation of pressures higher than CJ. It should be noted that such collisions of strong waves

would probably not be observed very often experimentally, since the symmetry imposed in the

computation would not occur. The effect is enhanced because the calculation is carried out in

two dimensions.

Further Discussion of the Computational Results

All of the behavior described above can be seen from a different perspective in Figures

15-18, which show the computed pressures, temperatures, turbulent kinetic energy, and rate

of consumption of hydrogen in a series of three-dimensional plots over 7 of the -chambers at5 different times. These plots illustrate the progress of the flame as it propagates down the

channel from an axial position of about 16.5 m to 22 m. The flame position is defined by the

temperature plots in Figure 16. Note that burning continues for several chambers behind the

flame front. Note also in Figure 17 that the flow ahead of the flame has generated turbulence,

35
primarily because of the strong shear in the flow behind the obstacles. Recall that the threshold

turbulent kinetic energy km is 4.74 J/kg. The turbulence levels shown in Figure 17 (peaks of

about 10 kJ/kg) far exceed the threshold level. This provides justification for introducing the

threshold as a startup mechanism to simply provide a minimum burn velocity.

Because of the choking effect, the flame does not actually accelerate rapidly in the region

shown in these figures until DDT occurs. However, the pressures are still rising rapidly in the

choked regime.. It is important to note the mechanism leading to this pressure rise. As the

pressure wave propagates down the channel, it hits the obstacles. Reflection of these shocks

off the front faces of the obstacles results in a large increase in pressure and some shock heating

(see Figures. 15 and 16). In accordance with Equation (27), this increases the burn rate (see

Figure 18), which feeds back into an increase in the overall pressure level of the wave.

There is a further mechanism contributing to this feedback in the computation. When the

gas burns in the corners in front of the obstacles, the increase in pressure there forces the gas

out past the corner of the obstacle and into the central flow. (See the flow pattern in Figure

19.) The resulting shear stresses contribute to a strong increase in the production of turbulence.

This, in turn, results in a further increase in the burn rate.

The computation exhibits the deflagration to detonation transition at about 315 ins, when

the flame has reached an axial position of 22 m. Note that a sharp pressure ridge has formed

all the way across the channel in the uppermost plot in Figure 15; i.e., the pressure spike is not

restricted to the region just in front of an obstacle.

As noted above, detonation was not observed experimentally at the concentration of 14.5%

H2 used in this calculation, but was observed at 15% H2. It is likely that the actual concentrations at which DDT is or is not experimentally observed would fluctuate from experiment to

experiment, so the computational prediction of a transition at 14.5% H2 for the particular m odeling parameters used is not regarded as being very conclusive. The fact is that we do not claim

to be able to simulate the detailed physical processes which occur in DDT. But it is reasonable to assume that the computed trends in pressure, temperature, and turbulence levels are

indicative of real processes, and that what is being calculated is the evolution of a system that

is progressing toward a tendency to detonate. (Note that it was assumed in Section III that the

pressure was constant through the flame when determining the dependence of the combustion

rate on temperature. This is not true for a detonation; however, the propagation velocity for a

detonation is determined by the conservation laws and driving forces, and not the combustion

rate, so this discrepancy is irrelevant for the numerical detonation, regardless of its validity on

other grounds.)

Results of Variations in the Computational Parameters

In order to assess the accuracy of the calculation, a second calculation was performed on

a refined grid. The baseline computation required about 31 hoursofcmuetienaCry

15S computer. Because of the expense of these calculations, it was not possible to reduce the

grid size enough to make a thorough study of convergence. What was done wa's to decrease

the spacing by a factor of 2/3. (The grid. for the complete domain was changed from 13 x 397

to 19 x 595.)

36
t

315 ms

312 ms

310 ms

307 ms

303 ms

2561

5.0'

0

90

x

Figure 15. Three-dimensional plots of pressure from the baseline calculation shown in Figures

10-14. Only the region between the 7th through the 14th obstacles along the FLAME facility

channel are shown. Coordinates are indicated in centimeters, and absolute pressure is given in

atmospheres. The configuration has a symmetry plane at x = 0. The times corresponding to

each plot appear on the right side of the figure. The axial positions of the pressure transducers

which provided the pressures shown in Figures 11-14 are indicated by bullets (.). (Note: In all

three-dimensional plots, the values of the dependent variables (pressure, temperature, turbulent

kinetic energy, and rate of combustion of hydrogen) have been set equal to zero in the region

occupied by the obstacles. -this is done simply to make the obstacles easily identifiable; the

variables are actually undefined there.)

37
t

315 ms

312 ms

310 ms

307 ms

303 ms 2561

T 2000

0

90

x

Figure 16. Three-dimensional. plots of temperature (in degrees Kelvin) from the baseline calculation.

38
to

315 ms

312 ms

310 MS

307 ms

303 ms

2561

10

k 5

0

90

x

Figure 17. Three-dimensional plots of turbulent kinetic energy (in kJ/kg) from the baseline

calculation.

39
315 ms

312 ms

310 ms

307 ms

2561 303, ms

Fuel rate

50

25

0

x

Figure 18. Three-dimensional plots of the hydrogen combustion rate (in kg/(M3-S)) 'from the

baseline calculation.

40
Vmax = 580 m/s

.1

I,...

I,..

.. tl

ii It.

.,ii I,,. .. tl It.. * I t* It

It...

It.

I,.. * I,.

I.... * .. ~ I.

I...

* ;r~ \\\ ,*

* i/Il/iT * - .- ////It \~~'..v..... * - ~ ~ -

* ~ I l~t\\..~

* ~----------

* - Ill ~ -

at. ~

ax ...... ''Ifliits.. --- 'a

tax ss.--~' ,~lills.. a a

a ax s....-' JIIIIIS.... a

S...-'

a a xxs.... .ai lilt ~

a ax xSS' ,tlltt~. a

it Sx.S' t lilt. -- -till, a a a *xxSS. - t tilt, *--'~ I I) a

a a a , 5 5 5 .. ht iiit.. ---- i, a a

a a' .... %x-.' ~

* a.. .. x~..- - titi,..-

-'a',-- -

ta..*'t~*t it.

-~------10th obstacle

(z = 18.3 m)

Figure 19. Computed flow pattern between the 9th and 11 th obstacles at t = 310 ins. The

calculation is actually restricted to the domain x'> 0, with a symmetry boundary condition at

x = 0. However, in this figure the velocity vectors are reflected across the midplane to show

the full width of the channel for clarity.

41
Some of the results are presented in Figures 20-22. The flame trajectory (see Figure 20)

is seen to be quite similar to that obtained from the previous calculation, although the flame

accelerates faster on the finer grid and detonates sooner. These differences are the result of

relatively small errors in the numerical solution on the coarse grid compared to the fine grid.

Because of the feedback, mechanisms described above, any errors in the flow variables or

.turbulence parameters become amplified.

To illustrate this sensitivity, an additional simulation was carried out (on the coarse grid)

in which the only difference from the baseline calculation was that the rate of combustion during

ignition was increased by 31%. The results are given by the dashed line in Figure 20. It is

seen that, even though ignition ends before the flame has propagated very far, the gas flow has

been given a boost sufficient to generate a higher initial turbulence level and cause considerably

greater acceleration of the flame. This example shows how susceptible the progress of such

experiments can be to initial conditions. If the burn in the first chamber proceeds differently

from shot to shot because of, e.g., different manifestations of fluid instabilities, the end results

may vary significantly. For example, if the concentration is marginally detonable, the startup of

the burn may determine whether detonation actually occurs.

The pressure histories in Figure 21 and the three-dimensional pressure-plot i'n Figure 22

display the behavior suggested earlier; when numerical dissipation is diminished by lowering

,the grid size, the shocks and, acoustic waves are resolved better and exhibit more structure.

Two more calculations were made to determine the sensitivity of the results to the parameters used in the combustion model. In the first, the values of Ac,km, and Lm were adjusted

so that the slope dST/dVik_ was reduced to 0.82. (See Figure 5.) The initial burn velocity STi

and the flame thickness b, were changed to 3.6 m/s and 17 cm, respectively. The new values

Of Ac,km, and Lm, were 18, 15.6 J/kg and 52 cm. For this run, the pressure and temperature

exponents m and n were held at 0.2 and 1.64, as before.

For the second calculation, these same new values of Ac,km,, and Lm were used, but the

temperature exponent n was raised to 2.0. The pressure exponent m remained at 0.2. Both

these calculations were made on the 13 x 397 grid.

The results are compared with those from the original computation in Figures 23-25.

Qualitatively, the effects on the flame trajectories are entirely predictable. The flame does not

burn as rapidly as before with either of the new sets of parameters, but a more r apid burn is

obtained for the higher value of n. In both cases, a numerical detonation occurs.

As seen in the sample pressure histories in Figure s 24 and 25, the pressure levels are

lower with the new parameters, also as expected.

V. CONCLUDING REMARKS.

This paper has presented a method for simulating large-scale premixed turbulent flames.

The emphasis here has been on developing a way to model the burn rate so that the gas

flow driven by the expansion of the burned gases can be computed. A crucial aspect of the

.application of the model to the computation of the behavior of accelerating flames is that the

resulting turbulence intensity and compression of the gas feeds back into the, combustion rate..

42
.40-

0

0.

CL

E

U..

30-

20-

10-

--. V-- Experiment F-2 3

.. .... .B a s e lin e c a ic ,. (coarse grid)

Fine grid

-----Fast ignition

S.0

a

7?

I'.

''a

I,' 'a'

~IFeI

U

.1

0 -b

0 100 200 300 400

Time (ins)

Figure 20. Comparison of experimental flame trajectories with those obtained from the baseline

calculation, from a similar calculation on a refined finite-difference grid, and from a calculation

in which the ignition strength was increased by 31 % (see text).

43
8 z = 18.6 m

Experiment (PKU5)

Experiment (PKU1O) ....... .. Caic. (fine grid) 6

0.

cc 4

a

2

0

290 295 300 305 - '310 315

Time (ins)

Figure 21. P~ressure histories at 18.6 m: From the experiment and from the calculation on the

refined grid.

44
305 ms

303 mS

301 ms

299 ms

297 ms

2561

P5.0

0

90

x

Figure 22. Three-dimensional plots of pressure from the calculation on the fine grid. The bullet

indicates the location of the pressure transducer corresponding to the results shown 'in Figure

21.

45
40- -7V-- Experiment F-2302

Baseline caic.

(S'T - 1.14, n = 1.64

C

0

cc

0

E

Lu

30-

20-

10-

.....S'T = 0.82, n = 1.64

.S)'T = 0.82, n = 2.0 a

77

is I,,

0

0 100 200 300 400

Time (ins)

Figure 213. Comparison of the flame trajectory obtained from the baseline calculation."and from

the calculations in whi ch the combustion model was altered (see text).

146
z = 18.6 m

8

Experiment (PKU5)

------ Experiment (PKU1lO) ......... Calculation (S'T 0.82., n .= ,;1.64)

6

EU

I I s

315 320 325 330 335- 340

Time (ins)

Figure 24. Pressure histories at 18.6 m: From the experiment and from the. calculation with

reduced combustion rate.

4.7
z = 18.6 m 8

6

Experiment (PKU5) ....... - Experiment (PKU1O)

............ Calculation (S'T = 0.82, n =2.0)

0-%

E

CL

£

S

4

2

01-

305 310 .315 320 325 330

Time (ins)

Figure 25. Pressure histories at 18.6 m: From the experiment and from the calculation with

reduced combustion rate, but with temperature exponent n = 2.

48
A second important aspect, pertinent to the simulation of large-scale flames, is that the flame

thickness is artificially increased and controlled for satisfactory numerical behavior.

The development of the model has been guided by the desire to use as much general

experimental information as-possible to determine parameters in the model, while minimizing

,the amount of knob-twisting to tune to specific experiments. The application to the sim'ulation.

of experiments in the FLAME facili ty has been quite successful. The simulated flame trajectories

and pressures provide reasonably accurate reproductions of experi mental behavior. -Some of the

computed pressure histories exhibit amplitudes and low frequencies which agree quantitatively

with experimental data to within about 25%. It has been sh .own that computation of pres .sures

which match experiment follows from a good reproduction of the rate of burnout of the individual

chambers formed by the obstacles. These rates, in turn, depend on the -accuracy of the burn

velocities. We have been able to obtain adequate burn velocities by requiring the .nurnerical

flames to behave in a way co .nsistent with known properties of turbulent flames and applying

flame speed data from the weakly turbulent regime of the experiments.

The resulting combustion model should be of interest to those confronted with the type

of simulation addressed here; viz., one in which the size of the computational domain greatly

exceeds the flame thickness. In order to apply the method, the following are required: (1) Some

experimental and/or theoretical information on turbulent flame speeds and their dependence on

the state of the gas, and (2) A computer code which solves the equations of reacting flow (for

a one-step reaction), computes the state of turbulence in the gas, and provides for flexibility in

the specification of the burn rate and diffusion coefficients. Once a numerical flame with the

desired properties has been developed, it can be used for the prediction of flame behavior in a

variety of situations.

These computations are not in tended to describe deflagration-to-detonation' transitions.

However, they do graphically demonstrate the way in which feedback mechanisms produce a

combustion wave which propagates down the FLAME channel with increasing gas velocities

and pressure levels. This suggests that there is a potential for the eventual development of the

ability to numerically simulate the evolution of the state of the gas as it approaches a sensitivity to

detonation. Given a suitable experimental data base, this might provide a means for predicting

DDT.

The combustion model should still be regarded as being under development. In particular,

there are two features specific to flame acceleration which need improvement. The first is that

a mechanism should be provided that will simulate the sudden burst of acceleration observed

in the experiments (the sharp knee in the experimental curve in Figures 10, 20, and 23). The

second is the tendency of the numerical flame to exhibit a transition to detonation in the absence

of corresponding experimental behavior. As noted earlier, however, this latter problem may not

be particularly severe, in view of the demonstrated detonability of mixtures at concentrations

close to that studied here.

One way to inhibit numerical DOT would be to include a mechanism for flame quench

due to turbulent flame stretch, which is a plausible physical process. Simulating the rapid

acceleration while still inhibiting detonation appears less straightforward. It may be that a more

accurate treatment of thermal mixing due to turbulence is required, perhaps in concert with

compressional heating, but this is merely speculation.

49
Future work should address these questions. In particular, a joint experimental and computation al study of the effects of varying hydrogen concentration would be of interest. Also of

interest to the reactor safety community would be the inclusion of gas'venting in the simulations.

It should. be emphasized that the mechanism for the extremely rapid acceleration of the flame

fronts which appears in the present work arises to a large extent from the fact that the channel

is enclosed except for the opening at one end. Pressure relief due to venting can reduce flame

speeds dramatically3'

5. However, in some cases, vent ing may enhance. the combustion rate,

presumably by introducing additional turbulence.

Another direction that future research could take would be to break the computational

problem with its large spatial domain up into smaller pieces. Two or three chambers could be

considered with appropriately chosen inlet and outlet boundary conditions. This would permit grid refinement to ascertain the convergence properties of the solution. It would also allow better

resolution of the turbulent flame fronts, recirculation regions and shock and acoustic waves.

50
APPENDIX A. INTEGRATION OF PRODUCTION-DISSIPATION TERMS

The following numerical scheme is used to ensure that numerical errors in the production

,and dissipation terms do not result in. nonnegative k or c. First, the production and dissipation

terms are split off from the rest of the two equations. This means that in each cell, these terms

are allowed to advance k and c within one computational time step according to

d-(pk) P - P dt

d C 2

-(P0) ci= E dt k

(A-i1)

(A - 2)

where

P = U:Vu

(Mesh subscripts are dropped for simplicity. Ordinary derivatives are indicated, since the

quantities are decoupled from those in other computational cells. The convection and diffusion

terms. in Equations (10) and (11) are added separately.)

In the strong turbulence regime we have

A'T > A'L

This is essentially always true in our calculations. Then

A C Pk2

Define

Q pk

E PE

Then

CAQ2

and

51
P =CIUfP U(U~,

where

f, (u,k,f) a:VU

The rationale behind separating out the term fp(u,k,c) is that it will be held constant while

equations (A-i) and (A-2) are integrated over one time step. This is justifiable in that f,, is not a

strong function of k and c. (In incompressible flow, it has no k or Edependence at all.)

Equations (A-i) and (A-2) can now be written

dQ _ __ E (A- 3)

dt -E

dE E2

__ = Ci-PC2,- (A - 4) dt Q

where

P' CMfp(u,k,c)

From (A-3) and (A-4) we obtain

d (Q) 1idQ Q dE.

dt E E dt E2 dt

--( Ci - I + (C2, 1)

Let QIE. Then

-ý (c1, - 1)PIC + (C2 ,

and we use the following implicit scheme to integrate ~

-.1 + (C 2 , - i)At-5

1+(C1, - 1)pt~nAt ( 5

where superscripts n and n + 1 denote old and new time levels and At is the time step. Since

ClE,C2, > 1, positive ý' insures positive en+1*

52
From Equation (A-3),

dQ Q -dand the scheme

Qfl+l = ( 1Ptn+)Q . ,(A.- 6)

is consistent and results in nonnegative Qn+'. Applying first (A-5), then (A-6), and' using

En~l- Qf+l

permits evaluation. of

k n+l Qfl+l/p

n+1~ Ef+l /p

and these quantities are positive as lo ng as kn and En are positive.

53
APPENDIX B. DIFFERENCE FORMULA FOR f AT WALLS

A function which varies as 1l/y cannot be differentiated accurately by the usual secondorder central difference. formulas. Since E has this behavior near walls, a special differentiation

formula was derived for this purpose, as described below. A formula which is exact for 1 /Y is

af 3 Ili+ Ei-Ei2

a2 4, Ay

where i is a grid index. (In Conchas-Spray, c is a cell-centered quantity, so i

9 iday is centered

on the cell faces.) We have found that such an alteration at only the first grid point in from the

wall yields a more accurate solution; however, it is an inconsistent formula for the' derivative,

and should not be carried into the interior of the fluid. The disc ontinuity in transferring from this

formula to the second-order approximation

9y ~ Ay

obviously leads to some inaccuracy, but it appears to be less severe then. carrying the secondorder approximation all the way to the wall.

Note added in proof: In the most recent version of the Kiva code (a descendant of

Conchas-Spray 17) , Los Alamos researchers are solving the C-equation by advecting the length

scale L rather than c. This alleviates the aforementioned difficulties with c at wall boundaries.

54
REFERENCES

11. Fisk, J. W., Sherman, M. P., Tieszen, S. R. and Benedick, W. B. "Experimental Facilities to

Study Hydrogen Combustion and Detonation." Designing for Hydrogen in Nuclear Power

Plants, Niyogi, K. K. and Bernstein, M. D., Eds. pp. 25-30.. A.S.M.E., New York, Y.Y.

(1984)

2. Sherman, M. P., "Hydrogen Flame Acceleration and Detonation." Twelfth Water Reactor

Safety Research Information Meeting, Oct. 1984. pp. 338-356, NUREG/CP-0058 Vol. 3,

U. S. Nuclear Regulatory Commission (1985).

3. Sherman, M. P., Tieszen, S. R., Benedick, W. B., Berman M., and Carcassi, M. "The

.FLAME Facility-The Effect of Transverse Venting on Flame Acceleration and Transition to

Detonation," in "Dynamics of Explosions," Progress in Aeronautics and Astronautics, vol.

106, J. R. Bowen, et al, eds., pp. 66-89 (1986).

4. Lee, J.H.S., Knystautas, R. and Freiman, A., "High Speed Turbulent Deflagrations and

Transition to Detonation in H2-Air Mixtures." Combustion and Flame 56, 227 (1984)..

5. Chan, C., Moen, 1. 0., and Lee, J. H. S., ",Influence of Confinement on Flame Acceleration

Due to Repeated Obstacles." Combustion and Flame 49, 27-39 (1983).

6. Lee, J. H., Knystautas, R., and Chan, C. K., "Turbulent Flame Propagation in ObstacleFilled Tubes," Proc. Twentieth Symposium (international) on Combustion, The Combustion Institute, Pittsburgh (1984).

7. Moen, 1. 0., Donato, M., Knystautas, R., and Lee, J. H., "Flame Acceleration Due to

Turbulence Produced by Obstacles," Combustion and Flame 39, 21-32 (1980).

8. Henrie, J. 0. and Postma, A. K. "Analysis of the Three Mile Island (TMI-2) Hydrogen Burn,"

in "Thermal-Hydraulics of Nuclear Reactors," Volume 11, Amer. NucI. Soc. LaGrange Park,

Ill., U.S.A.S.M.E. 1157-1170 (1983).

9. Jones, W. P., "Models for Turbulent Flows with Variable Density and Combustion," in "Prediction Methods for Turbulent Flows," W. Kollmann, ed., McGraw- Hill, New York (1980).

10. Rodi, W., "Turbulence Models and Their Application in Hydraulics," International Association for Hydraulic Research, Rotterdamseweg 185-P.O. Box 177, 2600 MH Delft, The

Netherlands (1980).

11. Magnussen, B. and Hjertager, B., "On Mathematical Modeling of Turbulent Combustion

with Special Emphasis on Soot Formation and Combustion," Proc. Sixteenth Symposium

(International) on Combustion, The Combustion Institute, Pittsburgh, pp. 719-729 (1976).

12. Spalding, D. B., "Mixing and Chemical Reaction in Steady Confined Turbulent Flames,"

Proc. Thirteenth Symposium (international) on Combustion, The Combustion Institute,

Pittsburgh, pp. 649-657 (1971).

13. Bray, K. C. N., "Turbulent Flows with Premixed Reactants," Chapter 4 of "Turbulent Reacting Flows," P. A. Libby and F. A. Williams, Eds., Springer-Verlag, Berlin (1980)

14. M~arx, K. D., Lee, J. H. S., and Cummings, J. C., "Modeling of Flame Acceleration in Tubes

with Obstacles," Proceedings of the 1 1th World Congress of the International Association

for Mathematics and Computers in Simulation, Oslo, Norway, August 5-9, 1985.

55
15. Thompson, J. F., (Ed.), "Numerical Grid Generation," North-Holland, Amsterdam (1982).

16. Smooke, M. D. and Kosykowski, M. L., "Two-Dimensional Fully Adaptive Solutions of

Solid-Solid Alloying Reactions," J. Comp. Phys. 62, 1 (1986).

17. Cloutman, L. D., Dukowicz, J. K., Ramshaw, J. D. and Amsden, A. A., "Conchas-Spray:

A Computer Code for Reactive, Flows with Fuel Spray," Los Alamos National Laboratory

Report LA-9294-MS, May, 1982.

18. Libby, P. A. and Williams, F.A. in"Turbulent Reacting Flows," loc. cit., Chapter 1.

19. Williams, F. A., "Combustion Theory," 2nd Edition, Addison-Wesley, Reading, MA ý(1985).

20. Harlow, F. H. and Amsden, A. A., "A Numerical Fluid Dynamics Calculation Method for

All Flow Speeds," J. Comp. Phys, 8, 197 (1971).

21. Tennekes, H., and Lumley, J. L.,"A first Course in Turbulence," MIT Press, Massachusetts

Institute of Technology, Cambridge, MA (1972).

22. Hinze, J. 0., "Turbulence," McGraw-Hill, New York (1975).

23. Butler, T. D., and O'Rourke, P. J., in Proc. Sixteenth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, pp. 1503-1516 (1976).

24. O'Rourke, P. J. and Bracco, F. V., "Two Scaling Transformations for the Numerical Computation of Multidimensional Unsteady Laminar Flames," J. Comp. Phys. 33, 185 (1978).

25. Marx, K. D., unpublished work on FLAME experiments without obstacles.

26. Libby, P. A., Bray, K. N. C., and Moss, J. B., "Effects of Finite Reaction Rate and Molecular

Transport in Premixed Turbulent Combustion," Combustion and Flame 34, 285 (1979).

27. A. N. Kolmogoroff, 1. G. Petrovsky, and N. S. Piscounoff, "Etudes de l'6quations de Ia

diffusion avec croissance de Ia quantit6 de mati6re et son application a un probl6me

biologique," Bull. Univ. Moscou. Ser. Internat., Sec. A, 1 No. 6, 1 (1937).

28. Clavin, P., and Li-nan, A., "Theory of Gaseous Combustion," NATO Advanced Study Institute Series B, Physics 116, Plenum Press, pp. 291-338 (1984).

29. Kerstein, A. R., "Pair-Exchange Model of Premixed Flame Propagation," to appear in Proc.

21st Symposium (International) on Combustion, Munich, Federal Republic of Germany,

August 3-8, (1986).

30. Warnatz, J., "Concentration-, Pressure-, and Temperature- Dependence of the Flame Velocity in Hydrogen -Oxygen -Nitrogen Mixtures," Comb. Sci. Tech. 26, 203 (1981).

31. Kee, R. J., Miller, J. A., and Jefferson, T. H., "CHEMKIN: A General-Purpose, ProblemIndependent, Transportable, Fortran Chemical Kinetics Code Package," Sandia. Laboratories Report SAND80-8003, Sandia National Laboratories, Livermore, CA (1980).

32. Hjertager, B., "Numerical Simulation of Turbulent Flame and Pressure Development in

Gas Explosions," in "Fluel-Air Explosions," University of Waterloo Press, p. 407 (1981).

33. Sherman, M. P., Tieszen, S. R., and Benedick, W. B. "The Effect of Obstacles on Flame

Acceleration and Transition to Detonation in a Large Channel," to appear in Proc. 21st

56
Symposium (international) on Combustion, Munich, Federal Republic of Germany, August

3-8, 1986.

34. Sherman, M. P., Tieszen, S. R., and Benedick, W. B., private communication (1986)

57
UNLIMITED RELEASE:

U. S. Government Printin g Office

Receiving Branch (Attn: NRC Stock)

8610 Cherry Lane

Laurel, MD 20707

(250 Copies for R3)

U. S. Nuclear Regulatory Commission (12)

Division of Accident Evaluation

Office of Nuclear Regulatory Research

Washington, DC 20555

Attn: B. W. Morris T. Lee

P. Worthington 3. Mitchell

B. Burson R. Meyer

W. S. Farmer J. Telford

M. Fleishman R. W. Wright

M. Silberberg C. W. Nilsen

U. S. Nuclear Regulatory Commission (8)

Office of Nuclear Reactor Regulation

Washington, DC 20555

Attn: V. Benaroya K. I. Parczewski

W. R. Butler Z. Rosztoczy

J. T. Larkins T. M. Su

R. Palla C. G. Tinkler

U. S. Department of Energy

R. W. Barber

Office of Nuclear Safety Coordination

Washington, DC 20545

U. S. Department of Energy (2)

Albuquerque Operations Office

Attn: J. R. Roeder, Director

Transportation Safeguards

J. A. Morley, Director

Energy Research Technology

For: C. B. Quinn

R. N. Holton

P.O. Box 5400

Albuquerque, NM 87185

D. M. Austin

office of Energy Research

U. S. Dept. of Energy

ER-7, MS G-226 Germantown

Washington DC 20545

Mr. J. Coleman

Office of Energy Research

U.S. Dept. of Energy

MS G-226 Germantown

Washington, DC 20545

Marvin Gunn, Jr., Manager

Energy Conversion Technology Proj.

Div. of Energy Conversion and

Utilization Technologies C-12

U. S. Department of Energy

1000 Independence Avenue

Washington, D.C. 20585

Ramesh Jain

Buildings and Community Systems

U. S. Department of Energy

1000 Independence Avenue

Washington, D.C. 20585

Dr. A. H. Laufer

Fundamental Interactions Branch

Division of Chemical Sciences

Office of Basic Energy Sciences

U. S. Department of Energy

Washington, D. C. 20545

Dr. Danny C. Lim

Industrial Programts

U. S. Department of Energy

1000 Independence Avenue

Washington, D.C. 20585

Dr. 0. Manley

Division of Engineering

and Geosciences

Office of Basic Energy Sciences

Department of Energy

Washington, D.C. 20545

Dr. R. S. Marianelli

Division of Chemical Sciences

U. S. Department of Energy

Washington, D.C. 20585

K. D. Smith, Program Coordinator

Energy Technology Division

DOE Albuquerque Operations Office

Albuquerque, NM 87115

58
Donald K. Stevens, (Acting)

Associate Director

Basic Energy Sciences

U. S. Department of Energy

Washington, D. C. 20545

Dr. F. Dee Stevenson

Office of Energy Research

U. S. Dept. of Energy

R7142, MS G-226 Germantown

Washington DC 20545

Argonne National Laboratory (2)

Attn: Dae Cho

R. Anderson

9700 South Cass Avenue

Argonne, IL 60439

Brookhaven National Laboratory (2)

T. Ginsberg

G.. Greene

Upton, NY 11973

Los Alamos National Laboratory (7)

Attn: F. J. Edeskuty

R. Gido

J. Carson Mark

G. Schott

M4. Stevenson

J. Travis

K. D. Williamson, Jr.

P.O. Box.1663

Los Alamos, NM 87545

A. A. Amsden

T-3, MS-B216

Los Alamos National Laboratory

Los Alamos, NM 87545

T. D,. Butler

T-3, MS-B216

Los Alamos National Laboratory

Los Alamos, NM 87544

J. K. Dukowicz

T-3, MS-B216

Los Alamos National Laboratory

Los Alamos, NM 87545

J. Ramshaw

T-3, MS-B216

Los Alamos National Laboratory

Los Alamos, NM 87545

Dr. Elaine Oran

Naval Research Laboratory

Code 4040

Washington, D.C., 20375

Oak Ridge National Laboratory (2)

Attn: A. P. Malinauskas

T. Kress

NRC Programs

P.O. Box X, Bldg. 4500S

Oak Ridge, TN 37831

Prof. Philip J. Smith

Brigham Young University

Department of Chemical Engineering

350 C.B.

Provo, UT 84602

Prof. A. Chorin

Mathematics Department

University of California

Berkeley, CA 94720

Prof. John W.

University of

Dep artment of

Berkeley, CA

Daily,

California, Berkeley

Mechanical Engineering

94720

Prof. A. K. Oppenheim.

Dept. of Mechanical Engineering

University of California

Berkeley, CA 94720

Robert F. Sawyer

Mechanical Engineering

University of California, Berkeley

Berkeley, CA 94720

Prof. W. A. Sirignano

Office of the Dean

School of Engineering

University of California, Irvine

Irvine, CA 92717

59
Prof. W. G. Kollmann

Mechanical Engineering

University of California

Davis, CA 95616

Prof. C. K. Law

Dept. of Mechanical Engineering

University of California

Davis, CA 95616

University of California @ Los Angeles

Nuclear Energy Laboratory (2)

Attn: Prof. I. Catton

Prof. D. Okrent

405 Hilgard.Avenue

Los Angeles, CA 90024

Prof. P. Libby

Dept. of Applied Mechanics and

Engineering Science

University of California at

San Diego

La Jolla, CA 92093

Professor S. S. Penner

Energy Center - B-010

University of California,San Diego

La Jolla, CA 92093

Prof. T. G. Theofanous

Chemical and Nuclear Engineerng Dept.

University of California

Santa Barbara, CA 93106

Josette Bellan

Jet Propulsion Laboratory

California Institute Tech.

4800 Oak Grove Drive

Pasadena, CA 91109

Prof. P. E. Dimotakis

Graduate Aeronautical Lab 301-46

California Inst. of Technology

Pasadena, CA 91125

A. Leonard

Graduate Aeronautical Lab 301-46

California Inst. of Technology

Pasadena, CA 91125

Prof. K. N. C. Bray

University of Cambridge

Silver St.

Cambridge, CB3 9EW

England

Prof. Norman Chigier

Carnegie-Mellon Un iversity

Mechanical Engineering Department

Pittsburgh, PA 15213

Prof. Melvyn C. Branch

University of Colorado

Mechanical Engineering

Boulder, CO 80309

Universita Degli Studi Di Pisa

Attn: M. Carcassi

Dipartmento Di Costruzioni

Meccaniche E. Nucleari

Facolta Di Ingegneria

Via Diotisalvi 2

56100 Pisa

ITALY

Dr. V. P. Roan

College of Engineering

University of Florida

Gainesville, FL 32611

Universit~t Heidelberg

Attn: Juergen Warnatz

Heidelberg

FEDERAL REPUBLIC OF GERMANY

Professor R. Strehlow

Dept of Aeronautics

University of Illinois

Urbana-Champaign, IL 61801

Imperial College of Science and

Technology

Attn: Dr. A. D. Gosman

Dept. of Mechanical Engineering

Exhibition Road

London SW7 2BX

UNITED KINGDOM

60
Prof. James H. Whitelaw

Fluids Section

Dept of Mechanical Engineering

Imperial College of Science and Techology,

Exhibition Road

London SW7 2BX

Great Britain

Dr. A. C. McIntosh

Fuel and Energy Dept.

Leeds University,

Leeds LS2 9JT

ENGLAND

McGill University

Attn: Prof. John H. S. Lee

315 Querbes

Outremont, Quebeck H2V 3W1

CANADA

Simon Engineering Laboratory

Attn: Prof. W. B'. Hall

University of Manchester

M139PL

UNITED KINGDOM

Judson R. Baron

Dept. of Aeronautics & Astronautics

Massachusetts Inst. Technology

Room 33-217

Cambridge, MA 02139

Massachusetts Institute of Technology

Attn: N. C. Rasmussen

Nuclear Engineering Department.

Cambridge, MA 02139

Claus Borgnakki

Dept. of Mechanical Engineering

University of Michigan

550 East University

Ann Arbor, MI 48109

Dr. G. M. Faeth

University of Michigan

Ann Arbor, MI 48109-2140

University of Michigan

Attn: Prof. M. Sichel

Department of Aerospace Engineering

Ann Arbor, MI 47109

University of Michigan

Nuclear Engineering Department

Ann Arbor, MI 48104

Northwestern University

Attn: Prof. S. G. Bankoff

Chemical Engineering Department

Evanston, IL 60201

Pro~f. Bernard J. Matkowsky

Dept. of Engineering.Sciences

and Applied Mathematics

Northwestern University

Evanston, IL 60201

Prof. Frediano V. Bracco

Princeton University

Department of Mechanical and

Aerospace Engineering

The Engineering Quadrangle

Princeton University

Princeton, NJ 08544

Prof. Forman A. Williams

Dept. of Mechanical and

Aerospace Engineering

Princeton University

Princeton, NJ 08544

Prof. Colin R. Ferguson

Purdue University

Dept. of Mechanical Engineering

North Grant Street

West Lafayette, IN 47907

Rensselaer Polytechnic Institute

Attn: Dr. J. E. Shepherd

Troy, NY i2180-3590

Prof. C. T. Bowman

Stanford University

Mechanical Engineering, Bldg. 520

Stanford, CA 94305

Prof. William Reynolds

Mechanical Engineering Bldg.

Stanford University

Stanford, CA 94305

61
Institute f~r Kernenergetik und

Energiesysteme (2)

Attn: M. Buerger

H. Unger

University of Stuttgart

Stuttgart

FEDERAL REPUBLIC OF GERMANY

Prof. R. W. Bilger

Mechanical Engineering

University of Sydney

New South Wales,'2006

Australia

Technische Universit~t Miinchen

Attn: Dr. H. Karwat

8046 Garching

FEDERAL REPUBLIC OF GERMANY

Texas A

Nuclear

College

& M University

Engineering Dept.

Station, TX 77843

Prof. J. R. Bowen, Dean

College of Engineering

University of Washington

Seattle,.WA 98185

0

Prof. David T. Pratt

University of Washington

Mechanical Engineering FU-10

Seattle, WA 98195

University of Wisconsin

Attn: Prof. M. L. Corradini

Nuclear Engineering Department

1500 Johnson Drive

Madison, WI 53706

Prof. Mitchell D. Smooke

Dept. of Mechanical Engineering

Yale University

New Haven, CT 06520

AERE Harwell (2),

Attn: J. R. Matthews, TPD

Didcot

Oxfordshire OX11 ORA

UNITED KINGDOM

Atomic Energy Ltd. (2)

Attn: D. Liu

D. Wren

Whiteshell Nuclear Research

Establishment

Pinawa, Manitoba

CANADA

Atomic Energy Canada Ltd.

Attn: P. Fehrenbach

Chalk River, Ontario KOJ 1JO

CANADA

Ms. Kathryn L. Barnes

British Gas pic

Research & Development Division

Midlands Research Station

Wharf Lane, Solihull

West Midlands B91 2JW

Battelle Columbus Laboratory (2)

Attn: R. Denning

M. Leonard

505 Kind Avenue

Columbus, OH 43201

Battelle Institut E. V. (3)

Attn: Dr. Werner Geiger

Dr. Guenter Langer

Dr. manf red Schildknecht

Am Roemerhof 35

6000 Frankfurt am Main 90

FEDERAL REPUBLIC OF GERMANY

Battelle Pacific Northwest Laboratory (2)

Attn: M. Freshley

G. R. Bloom

P.O. Box 999

Richland, WA 99352

Belgonucleaire S. A.

Attn: H. Bairiot

Rue de Champ de Mars 25

B-1050 Brussels

BELGIUM

CEC

Director of Research

Science & Education

Attn: B. Tolley

Rue De La Loi 200

1049 Brussels

BELGIUM

62
Dr. Pierre Joulain

Chimie Phsicque deLa

Combustion.ERA16O au CNRS

Domaine du Def fend

Mignaloux Beauvoir

86800 St..Julien l'Ars

France

Dr. Paul Clavin

Dept. de Combustion

Centre St. Jerome

Marseille Cedex 13 13397

France

Paul Thibault

Combustion Dynamics Ltd.

Box 175

Ralston, Alberta

Canada TOJ 2N0

Dr. Rudy Maly

Daimler Benz AG

Abteilung E6T/E

Postfach 202

D-7000 Stuttgart

West Germany

Duke Power Company (2)

Attn: F. G. Hudson

A. L. Sudduth

P.O. Box 33189

Charlotte, NC 28242

John H. Pohl

Energy & Environmental Research Corp.

18 Mason

Irvine, CA 92718

EG&G Idaho (3)

Willow Creek Building, W-3

Attn: D. Croucher

R. Hobbins

Sever Sadik

P.O. Box 1625

Idaho Falls, ID 83415

Electric Power Research Institute (4)

Attn: J. Haugh

W. Loewenstein

B. R. Sehgal

R. Vogel.

3412 Hillview Avenue

Palo Alto, CA 94303

ENEA Nuclear Energ Alt Disp (2)

Attn: P. L. Ficara

G. Petrangeli

Via V. Brancati

00144 Roma

ITALY

ENSMA

Prof. Numa Manson

Rue Guillaume VII

86034 Poiti~ers,

FRANCE

Dr. Thomas S. Bustard

Energetics, Incorporated

9210 Route 108

Columbia, MD 21045

Factory Mutual Research Corporation

Attn: R. Zalosh

P.O. Box 688

Norwood, MA 02062

Fauske & Associates

Attn: R. Henry

16W070 West 83rd Street

Burr Ridge, IL 60521

James A. Kezerle

Gas Research Institute

8600 West Bryn Mawr Ave

Chicago, Il 60631

General Electric Corporation

Attn: K. W. Holtzclaw

175 Curtner Avenue

Mail Code N 1C157

San Jose, CA 95125

General Electric Corporation

Attn: M. I. Temme, Manager

Probabilistic Risk Assessment

Advanced Reactor Systems Dept.

P.O. Box 3508

Sunnyvale, CA 94088

Charles A. Amann

G. M. Research Laboratories

Twelve Mile and Mound Roads

P. 0. Box 9055

Warren, MI 48090-9055

63
Dr. Roger B. Krieger

G. M. Research Laboratories

Twelve Mile and Mound Roads

P. 0. Box 9055

Warren, MI 48090-9055

Gesellschaft fiur Reaktorsicherheit (GRS)

.Postfach 101650

Glockengasse 2

5000 Koeln 1

FEDERAL REPUBLIC OF GERMANY

Gesellschaft f~r Reaktorsicherheit (2)

Attn: Dr. E. F. Hicken

Dr. H. L. Jahn

8046 Garching

Forschungsgelande

FEDERAL REPUBLIC OF GERMANY

Institute of Nuclear Power Operation

Attn: Henry Piper

1100 Circle 75 Parkway, Suite 1500

Atlanta, GA 30339

International Technology Corporation

Attn: Mario H. Fontana

575 Oak Ridge Turnpike

Oak Ridge, TN 37830

Japan Atomic Energy Research Institute

Attn: Dr. K. Soda, Manager

Chemical Engineering Safety Laboratory

Dept. of Nuclear Fuel Safety

Tokai-mura, Naku-gun Ibaraki-ken

319-11

JAPAN

Japan Atomic Energy Research Institute

Attn: Dr. T. Fujishiro, Manager

Dept. of Fuel Safety Research

Tokai-mura, Naka-gun, Ibaraki-ken

319-11

JAPAN

Japan Atomic Energy Research Institute

Attn: Dr. Kazuo Sato, Director

Dept. of Reactor Safety Research

Tokai-mura, Naka-gun Ibaraki--ken

319-11

JAPAN

Kraftwerk Union

*Attn: Dr. M. Peehs

Hammerbacherstrasse 12 & 14

Postfach 3220

D-8520 Erlangen 2

FEDERAL REPUBLIC OF GERMANY

Lehrgebiet f~r Mechanik der

RWTH Aachen

Attn: Prof. Dr. Ing. N. Peters,

Templergraben 55

D5100 Aachen

FEDERAL REPUBLIC OF GERMANY

Dr. Ingar 0. Moen

Head/Shock and Blast Group

Defence Research Establishment

Suffield

Ralston, Alberta

Canada TOJ 2N0

Mississippi Power & Light

Attn: S. H. Hobbs

P.O. Box 1640

Jackson, MS 39205

National Nuclear Corp. .Ltd.

Attn: R. May

Cambridge Road

Whetestone, Leicester, LE8 3LH

UNITED KINGDOM

Netherlands Energy Research

Attn: K. J. Brinkmann

P.O. Box 1

1755ZG Pet ten NH

NETHERLANDS

NUS Corporation

Attn: R. Sherry

4 Research Place

Rockville, MD 20850

Dr. Arthur A. Boni

Physical Sciences Inc.

P. 0. Box 3100

Andover, MA 01810

Power Authority State of NY

10 Columbus Circle

New York, NY 10019

Foundation

64
Power Reactor Nuclear Fuel

Development Corp. (PNC)

Attn: Dr. Watanabe

FBR Project

9-13, 1 Chome, Akasaka

Minato-Ku, Tokyo

JAPAN

Projekt Schneller Brueter (3)

Attn: Dr. Kessler

Dr. Heusener

Dr. M4. Reimann

Kernforschungszentrum Karlsruhe

Postfach 3640

7500 Karlsruhe

FEDERAL REPUBLIC OF GERMANY

Pak-Yan Liang

Rockwell Int'l, Rocketdyne Div.

6633 Canoga Ave.,.FB47

Canoga Park, CA 91303

C. H. Priddin

Mgr. Combustion Methods

Rolls-Royce plc

P. 0. Box 31

Derby DE28BJ

Great Britain

Shell Research Ltd.

Paul H. Taylor

Thornton Research Centre

P.O.,Box 1

Chester CHl 3SH

ENGLAND

Jana Backovsky

SRI International

33 Ravenswood Avenue

Menlo Park, CA 94025

Stratton & Associates, Inc.

2 Acoma Lane

Los Alamos, NM 87544

Attn: W. Stratton

Statens Karnkraftinspektion

Attn: L. Hammar

P.O. Box 27106

S-10252 Stockholm

SWEDEN

Studsvik Energiteknik AB

Attn: K. Johansson

S-611 82 Nykoping

SWEDEN

.Swedish State Power Board

Attn: Wiktor Frid

S-162 Fach 87 Vallingby

SWEDEN

IJKAEA Safety & Reliability Directorate (2)

Attn: 3. H. Gittus

M. R. Hayns

Wigshaw Lane, Culcheth,

Warrington WA3 4NE

Cheshire

UNITED KINGDOM

UKAEA, Culham Laboratory

Attn: F. Briscoe

Abingdon

Oxfordshire 0X14 3DB

UNITED KINGDOM

UKAEA AEE Winfrith (6)

Attn: M. Bird, S. Board

T. Butland, D. Fletcher

R. Potter, A. Wickett

Dorchester

Dorset DT2 8DH

UNITED KINGDOM

Westinghouse Corporation (3)

Attn: N. Liparulo

J. Olboeft

V. Srinivas

P.O. Box 355

Pittsburgh, PA 15230

L. Cloutman, LLNL, L7035

C. K. Westbrook, LLNL, L-321

Gladys Shaw (15)

Sandia National Laboratories

P.O. Box 5800

Albuquerque, NM 87185

65
1131

1510

1512

1513

1530

1636

6000

6400

6412

6415

6419

6420

6422

6425

6425

6427

6427

6427

6427

6427

6427

6427

6427

6427

6427

6440

800.0

8100

8200

8240

8244

8245

8300

W. B. Benedick

J. W. Nunziato

J. C. Cummings

M. R. Baer

L. W. Davison

Attn: 1534 J. R. Asay

F. G. Blottner

D. L. Hartley

A. W. Snyder

A. L. Camp

F. E. Haskin

K. D. Bergeron

J, V. Walker

D. A. Powers

W. J. Camp

M. F. Young

M. Berman

D. F. Beck

J. T. Hitchcock

B. W. Marshall, Jr.

L. S. Nelson

M. P. Sherman

S. E. Slezak

D. W. Stamps

S. R. Tieszen

C. C. Wong

D. A. Dahigren

3. C. Crawford

E. E. Ives

R. J. Detry

C. W. Robinson

Attn: 8244 C. M. Hartwig

8245 R. J. Kee

8245 A. E. Lutz

D. R. Cheoweth

S. Paolucci

P. L. Mattern

.Attn 8310 R. W. Rohde

8340 W. Bauer

830J. S. Binkley

Attn: 8351 J. Y. Chen

8351 R. P. Lucht

8353 G. A. Fisk

8354 R. E. Palmer

8357 R. J. Carling

8360 W. J. McLean

8361 D. R. Hardesty

Attn: 8361 T. Fletcher

8362 T. M. Dyer

Attn: 8362 T. T. Bramlette

8363 B. R. Sanders

Attn: 8363 W. T. Ashurst

8363 P. K. Barr

8363 H. A. Dwyer

.8363 A. R. Kerstein

8363 S. B. Margolis

8363 K.,D. Marx (7)

8400 R. C. Wayne

8265 Publication Division/

'Tech. Library Proc. Div. 3141

8024 P. W. Dean,

Central Technical Files (3)

3141 Tech. Library Proc. Div.(3

3151 W. L. Garner

1!

66
NRC FORM 335 U.S. NUCLEAR REGULATORY COMIMISSION 1. REPORT NUMBER lAuigned bv 1`I0C. add V64 No.. if env)

1243 2. NUREG/CR -4855- 3201.3202 BIBLIOGRAPHIC DATA SHEET SAND87-8203

SEE INSTRUCTIONS ON THE REVERSE. R3

2..TITLE AND SUBTITLE 3. LEAVE BLANK

ýDEVELOPMENT AND APPLICATION OF A COMPUTER MODEL FOR

LARGE-SCALE FLAME ACCELERATION EXPERIMENTS 4. DATE REPORT COMPLETED

MONTH YEAR

S. AUTHORISI April 1987

6. DATE REPORT ISSUED

K. D. Marx .MONTH YEAR

_____________________________________July 1987 7. PERFORMING ORGANIZATION NAME AND MAILING ADDRESS (Include 2. Codwl S. PROJECTITASKIWORIC UNIT NUM BER

Sandia National Laboratories ________________

Livemore 45509. CA FIN OR GRANT NUMBER

Division 8 .363A14

10. SPONSORING ORGANIZATION NAME AND MAILING ADDRESS fin-dude Zop Code) I Is. TYPE OF REPORT

Division of Reactor Systems Safety Technical

Office of Nuclear Regulatory Research ________________

U. S. Nuclear Regulatory Commission .b. PERIOD COVERED fauOCJU'dataS)

Washington, DC 20555

12. SUPPLEMENTARY NOTES

13. ABSTRMACT. (200 word or Wall

A new computational model for large-scale premixed flames is deve loped and applied to

the simulation of flame acceleration experiments. The primary objective is to

circumvent the necessity for resolving turbulent flame fronts; this is imperative

because of the relatively coarse-computational grids which must be used in

engineering calculations. The essence of the model is to artificially thicken the

flame by increasing the appropriate diffusivities and decreasing the combustion rate,

but to do this in such a way that the burn velocity varies with pressure,

temperature, and turbulence intensity according to prespecified phenomenological

characteristics. The model is particularly aimed at implementation in computer codes

which simulate compressible flows. To this end, it is applied to the two-dimensional

simulation of hydrogen-air flame acceleration experiments in which the flame speeds

and gas flow velocities attain or exceed the speed of sound in the gas. It is shown

that many of the features of the flame trajectories and pressure histories in the

experiments are simulated quite well by the model. Using the comparison of

experimental and computational results as a guide, some insight is developed into the

processes which occur in such experiments.

14, DOCUMENT ANALYSIS - a.KEYWORDSIOESCRiPTORS

b. IDENTIFIERSIOPEN.ENDED TERMS

67'


Org. Bldg. Name Rec'd by jorg. Bldg. Name Rec'd by

Org. Bldg. Name Rec'd by jorg. Bldg. Name Rec'd by

I __

__ I __ __

~JSandia National Laboratories