^AT-L INST, OF STAND & TECH R.I.C.

Soft Furnishings by a Cigarette

U.S. Department of Commerce

Technologv Administration

National Institute of

Standards and TECHNOt-oev

August 1993

IMIST Special Publication 852

Nisr

-QG100

.U57

#852-5

1993
Fire Safe Cigarette Act of 1990

Under the Cigarette Safety Act of 1984 (P.L. 98-567), the Technical Study Group on Cigarette and Little

Cigar Fire Safety (TSG) found that it is technically feasible and may be commercially feasible to develop

a cigarette that will have a significantly reduced propensity to ignite furniture and mattresses. Furthermore, they found that the overall impact of such a cigarette on other aspects of the United States

society and economy may be minimal.

Recognizing that cigarette-ignited fires continue to be the leading cause of fire deaths in the United

States, the Fire Safe Cigarette Act of 1 990 (P.L. 101-352) was passed by the lOlst Congress and signed

into law on August 1 O, 1990. The Act deemed it appropriate for the U.S. Consumer Product Safety

Commission to complete the research recommended by the TSG and provide, by August 1 O, 1993. an

assessment of the practicality of a cigarette fire safety performance standard.

Three particular tasks were assigned to the National Institute of Standards and Technology's Building

and Fire Research Laboratory:

• develop a standard test method to determine cigarette ignition propensity,

• compile performance data for cigarettes using the standard test method, and

• conduct laboratory studies on and computer modeling of ignition physics to develop valid,

user-friendly predictive capability.

Three tasks were assigned to the Consumer Product Safety Commission:

• design and implement a study to collect baseline and follow-up data about the characteristics of

cigarettes, products ignited, and smokers involved in fires,

• develop information on societal costs of cigarette-ignited fires, and

• in consultation with the Secretary of Health and Human Services, develop information on changes

in the toxicity of smoke and resultant health effects from cigarette prototypes.

The Act also established a Technical Advisory Group to advise and work with the two agencies.

This report is one of six describing the research performed and the results obtained. Copies of

these reports may be obtained from the U.S. Consumer Product Safety Commission,

Washington, DC 20207.
ng the Ignition of

Soft Furnishings by a Cigarette

Henri E. IVIitler

George N. Wsl-ton

Building and Fire Research Laboratory

National Institute of Standards and Technology

Gaithersburg, MD 20899

U.S. Department of Commerce

Ronald H. Brown, Secretary

Technology Administration

Mary L. Good, Under Secretary for Technology

National Institute of Standards and Technology

Arati Prabakhar, IDirector

NIST Special Publication 852

August 1993
National Institute of Standards and U.S. Government Printing Office

Technology Washington: 1993

Special Publication 852

Natl. Inst. Stand. Technol.

Spec. Pub. 852

169 Pages (Aug. 1993)

CODEN: NSPUE2

For sale by the Superintendent of

Documents

U.S. Government Printing Office

Washington, DC 20402-9325
TABLE OF CONTENTS

Page

List of Tables vi

List of Figures vii

I. INTRODUCTION 1

n. SUBSTRAT, A MODEL OF A SUBSTRATE SUBJECTED TO A

MOVING HEAT SOURCE 3

A. Introduction 3

B. Ignition Dynamics 4

1. Surface Heat Source 4

2. Heat Transfer 6

3. Time to Ignition 8

4. Pyrolysis of the Substrate 8

5. Gas Diffusion 10

6. Boundary Conditions for Gases 11

7. Oxygen Diffusion 13

C. Numerics 17

1. Introduction 17

2. Boundary Conditions 20

3. Air Gap 21

4. Variable Grid 22

5. Variable Thermal Properties 23

6. Pyrolysis 23

7. Time Integration 24

D. Experiments 26

1. Experimental Arrangement 26

2. Material Data 29

E. Results and Discussion 34

1. Inert Substrate 34

2. Results witii Pyrolysis 35

F. Summary for Section n 47

m. CIGARET, A MODEL OF A QUIETLY SMOLDERING, ISOLATED

CIGARETTE 49

A. General Description 49

iii
B. Dynamics of a Smoldering Cigarette 49

1 . Pyrolysis Rates 49

2. Stoicliiometry , 51

3. Regression (Burning) Rate 52

4. Reaction Rates 53

5. Heat Production Rate 53

6. Energy Balance 53

7. Distributions Within the Cigarette 54

8. Pressure 55

9. Gas Transport 57

10. Diffusion Coefficients 57

1 1 . Conduction Losses 59

C. Previous Modeling Efforts 59

D. Modeling the Cigarette 62

1. Assumptions 62

2. Governing Equations 64

3. Choices for CIG25 70

E. Numerics 71

1. Discretization of the Equations 71

2. Method of Solution 72

F. Improvements in CIGARET Over CIG25 72

1. Physics 73

2. Input and Output 73

3. Documentation 74

G. Results 75

1. Sensitivity of Calculation 75

2. Velocity of Smolder Wave 75

3. Temperature Profiles 77

H. User's Guide 79

1. Running CIGARET 79

2. Input and Output Files 79

3. Restarting a Run 80

4. Contents of the Data File 81

5. Producing a Data File 82

I. Summary of Section III 86

iv
IV. SIMULATING A BURNING CIGARETTE ON AN IGNITABLE SUBSTRATE 87

A. Introduction 87

B. Qualitative Description 87

C. Use of the Two Programs 88

D. Detailed Calculations of Interactions 89

1. Conductive Flux to Substrate 89

2. Radiative Flux to Substrate 91

3. Complete Flux 92

4. Effects of the Substrate on the Cigarette 92

V. ACKNOWLEDGEMENTS 95

VI. REFERENCES 96

APPENDIX A: Conduction Algorithm Tests A-1

APPENDIX B: TMPSUB2 and SUBSTRAT Users' Guide B-1

APPENDIX C: Analysis of Numerical Errors Produced by

a Runaway Reaction Rate C-1

APPENDIX D: Analysis of Ignition Experiment D-1

APPENDIX E: Thermophysical Data for Cotton E-1

APPENDIX F: Heat Transfer Coefficient in the Presence

of a Substrate F-1

APPENDIX G: Calculation of the Convective View Factor G-1

APPENDIX H: Longitudinal Dependence of the Radiation Flux

to the Substrate H-1

APPENDIX I: Transient Depression of Cigarette Temperature

Upon Deposition on a Substrate I-l

APPENDIX J: Mass Transfer Coefficients J-1

APPENDIX K: Nomenclature K-1

V
LIST OF TABLES

Page

1. Features of the Substrate Model 4

2. Substrate Ignition Delay Times 29

3. Comparison of Substrate Ignition Delay Times 45

4. Surface Temperatures Which Would be Attained by the Substrate at the

Measured Ignition Times If the Substrate Were Inert 47

5. Ignition Temperatures at a 100 °C/s Rate of Rise 47

6. Position of Cigarette Smolder Front at 60 Seconds as a Function of the

Choice of Isotherm 77

A-1. Grid Spacing Tests A-7

A-2. Time Step Control Tests A-8

A-3. 3-D Transient Conduction Tests A-9

D-1. Measured and Inferred Quantities From the Ignition Experiments D-3

vi
LIST OF FIGURES

Page

1. One-dimensional conduction 17

2. Air gap heat transfer model 21

3. Substrate coordinate system 23

4. Schematic of heat source for ignition tests 27

5. Flux profile from heat source, as measured by a total heat flux gauge 28

6. Time required to ignite the fabric, for different (initial) heat

flux exposures 30

7. Peak surface temperatures of substrate as a function of time,

for the four exposures 36

8. Peak surface temperatures of substrate exposed to Q = 25 kW/m^,

as a function of time 37

9. Temperature and density of central cell, adjacent

cells, and cells in next ring around center, as functions of time 38

10. Temperature of central surface cell (i.e., peak temperature)

for three different grid sizes: 1.25, 0.50, and 0.25 mm cubes 39

11. Peak temperature for the 25 and 34.5 kW/m^ cases, assuming

(a) no pyrolysis, and (b) all three pyrolytic reactions 40

12. Peak temperature for the 18 kW/m^ case, with the same

assumptions as in Figure 11 42

13. Peak temperature for the 18 kW/m^ case, for several values of mean O2 mass fraction 43

14. Peak temperature for the 18 kW/m^ case, with various assumptions for

the pyrolysis 44

15. Peak temperature as a function of time, for all four cases, using

the best set of input data 46

16. Schematic of a smoldering cigarette 50

vii
17. Isotherms in the gaseous part of the cigarette during quiet, steady

burning 50

18. Oxygen concentration in the freely smoldering cigarette during quiet,

steady burning 56

19. Isotherms in the solid part of the cigarette during quiet, steady burning 56

20. Schematic of how air enters the smoldering cigarette when it is resting on a surface 58

21. Schematic of the energy losses of a smoldering cigarette as a function

of time, when it is dropped onto a surface 58

22. Result of a sample run with CIGARET 76

23. Longitudinal temperature distributions, along the axis and along the

surface, at t = 60 s, for a sample run with CIGARET 78

24. Flux emitted by a cigarette toward the substrate, along the contact line 90

B-1. TMPSUB2 Programs and Files B-1

B-2. Substrate Coordinate System and Variable Grid B-14

D-1. Total heat flux impinging on gauge, with and without purge flow D-2

D-2. Calculated heat transfer coefficient, as a function of the

temperature Tg of the impinging purge gas jet D-5

E-1. Thermal conductivity of cotton as a function of temperature

as measured by difl'erent workers E-4

F-1. Schematic of the exchange of energy between point P on the cigarette

surface and the substrate surface F-3

H-1. The geometric relationship between a general point on the cigarette

and another on the substrate H-2

viii
Certain commercial equipment, instruments, materials, or products are identified in this paper in order

to specify the experimental procedure adequately. Such identification is not intended to imply

recommendation or endorsement by the National Institute of Standards and Technology, nor is it intended

to imply that the materials or equipment identified are necessarily the best available for the purpose.

ix

MODELING THE IGNITION

OF SOFT FURNISHINGS

BY A CIGARETTE

Henri E. Mitler and George N. Walton

ABSTRACT

This paper describes the user-friendly computer models CIGARET and SUBSTRAT.

CIGARET calculates the time-dependent behavior of a cigarette smoldering quietly in the air, away from surfaces. The model incorporates diffusion and convection of gases, as well as the

kinetics of char oxidation. It calculates the internal heat fluxes, as well as the internal distribu- tions of temperature, gas velocity, and oxygen concentration. SUBSTRAT determines whether

a two-layer solid (with an air gap in between), exposed to a moving heating flux such as is produced by a cigarette, will ignite. Among the processes taken into consideration tire three- dimensional heat conduction in the substrate and its pyrolysis. This model has successfully simulated the thermal runaway signifying smoldering ignition of the substrfite when it is exposed to a set of external heating fluxes. SUBSTRAT and CIGARET have been designed

to work in tandem to simulate the most frequent cause of fatal flres: cigarette ignition of upholstered furniture and bedding. Users' guides are included.

Key words: cigarettes; cigarette model; computer model; free smolder; furniture fires; ignition; mathematical modeling; modeling; pyrolysis; simulation; smoldering; substrates

I. INTRODUCTION

Lighted tobacco products (cigars, cigarettes, pipes) continue to be the leading ignition source for fatal

fires in the United States (Miller, 1991). Most of diose fires are established in soft furnishings: uphol-stered furniture and bedding. As a result. Congress passed the Cigarette Safety Act of 1984, with a goal

to determine whether less fire-prone cigarettes were technically feasible. Under that Act, NIST developed

prototype mathematical models of the components of this ignition scenario. The reason for having this capability is to "test" the effects that some particular change(s) in a cigarette will have on its propensity

to ignite specified upholstered furniture quickly, cheaply, and repeatedly. The result was two prototype

programs, TEMPSUB (for TEMPerature of a SUBstrate) and CIG25 (Gann et al, 1988, Section 5).

Subsquent legislation, the Fire Safe Cigarette Act of 1990 (P.L. 101-352), mandated completion of such

research, namely to:

"(1) develop a standard test method to develop cigarette ignition propensity,

(2) compile performance data for cigarettes using the standard test method developed under

paragraph (1), and

1
(3) conduct laboratory studies on and computer modeling of ignition physics to develop valid,

"user-friendly" predictive capability."

This publication describes the research performed and the results obtained in responding to the third task. The completion of the first two tasks is described in "A Computer Model of the Smoldering Ignition of

Furniture," NISTIR 4973.

In this work, the ignition process is simulated using two computer codes, which are then used in tandem.

While a single code of a fully interactive cigarette and substrate is ideal, the physical and numerical

complexity of the combined burning near and at ignition renders it impossible to perform such

computations on the personal computers of today.

Section II describes SUBSTRAT, a model of the time-dependent heating of substrate when subjected to a moving heat source, and establishes a criterion for its ignition. As in most soft furnishings, the

substrate consists of a fabric layer over padding, with a potential air gap in between.

In order to examine how changing one or more properties of the cigarette will influence its ignition

propensity, it is necessary to understand its behavior when smoldering. This includes knowing how its external heat flux and burning velocity depend on its various geometrical, physical, and/or chemical

properties, and how these processes are modified when the cigarette lies on the substrate. Section III describes CIGARET, a model of a burning cigarette lying on a surface such as described above,

describable by SUBSTRAT. Both the physics and chemistry of tobacco pyrolysis and simultaneous heat

and gas transport in the cigarette are expressed by a set of coupled, nonlinear, partial difl'erential equations with nonlinear boundary conditions. The model produces the time-dependent distributions of

temperature, oxygen concentration, gas velocity, and burning rate in a freely smoldering cigarette which

has user-prescribed properties.

In Section IV, we then describe how to use the two programs to simulate a burning cigarette on a

susceptible substrate. SUBSTRAT reads an input file created by CIGARET, and also produces an output

file which is used by CIGARET.

2
n. SUBSTRAT, A MODEL OF A SUBSTRATE SUBJECTED TO A

MOVING HEAT SOURCE

A. Introduction

Consider what is involved in developing a computer model: we must simulate the behavior of a (typical)

substrate when it is subjected to a heat flux. In order to do that, the physics and chemistry of the process

must first be understood. It is then expressed as a set of equations describing the behavior; these

equations must then be solved, with the appropriate data and boundary conditions. The solution method

is numerical, and the data input is to a computer.

This program calculates the temperature of the upholstered furniture as a function of time and position,

when it is exposed to a prescribed local heating flux. The furniture is simulated as a flat, horizontal

cushioned seat; that is, a cushion consisting of fabric-covered foam padding. An air gap may be inserted

between the two layers. The prescribed flux can be given either by a simple Gaussian function, or via an input file; if it is the former, the flux can be highly peaked at a point, vary with time, and move at a constant (specified) rate over the top surface of the furniture, assumed to be horizontal. The radiative

and convective heat losses from the surface are given correctly. The temperature distribution within the

substrate is calculated, since the temperature history will determine whether ignition takes place. If and

when the temperature at a given location accelerates to a sufficiently high value (« 500 °C), we can say

that smoldering ignition has occurred. This is generally referred to as the ignition temperature, which

is not unique, even for a well-defined material, depending also on the magnitude of the incident heat flux. The ambient oxygen level can be set at whatever value one wishes.

There are some limitations to this program. It will not tell whether flaming ignition takes place. It also

does not treat the case where the flux is applied in a crevice, such as is formed between the seat cushion

and the seat back. The program does not take oxygen diflusion within the cushion explicitly into account; hence in certain threshold situations, where a small change in oxygen concentration determines whether

ignition does or does not take place, the results are ambiguous and not to be trusted. Note that it is often

difficult to obtain the needed kinetic and/or thermophysical parameters for the materials; or, when

available, to know how accurate they are. Therefore this caveat must also be made: even if the program

were perfect, its results are only as good as the input parameters which are supplied. On the other hand,

it should accurately reproduce or predict trends.

The first task is to examine the thermal response of the cushion to a specified flux distribution. This is done in Section II. B. [The flux data from the cigarette is described in Section III.] In Section II.C, the

numerical solution method is presented, and some of the physics is revisited in this (numerical) context. In Section II. D we discuss a set of experiments where a mock-up was ignited, give the measuroi ignition

times, and give the input data required to make computer runs. The results of these runs are given and

then discussed in Section II. E, validating the model to some degree. Section II.F summarizes the work.

For the reader interested in a more extensive discussion of this type of modeling, see the references cited

here, especially the bibliography in Gann et al. (1988). A copy of the program on disk may be obtained

from the authors, as well as a listing of the program and the program itself.

3
B. IGNITION DYNAMICS

The equations describing the relevant physical processes taking place in the substrate will be written in

this Section. They will mostly be written in general form to start, and then particularized to our problem.

The physical phenomena included in the model are summarized in Table 1 and discussed below.

Table 1. Features of the Substrate Model

Included Not Included

two porous layers with an air gap oxygen diffusion

3-d heat conduction melting and/or regression of the foam

variable (prescribed) thermophysical properties

correct (nonlinear) boundary conditions

endothermic (non-oxidative) pyrolysis

exothermic (oxidative) pyrolysis

char oxidation

arbitrary (prescribed) moving heat source

variable grid on a moving coordinate system

(approximate) radiative heat transfer within the material

impinging air flow, when wanted

1. Surface Heat Source

The potential ignition process is initiated by external heating of the surface. Thus, it is necessary to specify the heat gains at the bounding surfaces of the cushion. The cushion is assumed to be a rectangular parallelepiped, so a Cartesian coordinate system is employed. For the top surface, we first treat the localized heat flux from the glowing tip of the cigarette or other heat source. The flux consists

of two parts, one due to convection (0J and one due to radiation (0,):

<^u, is the heating flux reaching the surface. For a cigarette, these fluxes are strong functions of position,

since the glowing tip is only a few millimeters in extent. With this formulation, <t>j can be expressed in

the approximate form:

4
(2)

where Q = Q(r) is the sh^e, or view, factor of the cigarette as seen by the substrate at the point r. It

is only approximate in this form, because the cigarette surface temperature, Tj.jg, varies strongly with

position.

The convective flux is given by

cC^f.y.O = h[T^fg(x,yj) - r,(x.y,f)] (3)

where h is the heat transfer coefllicient, Tg is the surface temperature of the substrate, which is explicitly

shown to be a function of position and time, and where T^^^ is the cigarette surface temperature, also a

function of position and time. It is important to point out here that in the TMPSUB2 version of the

substrate model, the source flux must be specified by the user; in SUBSTRAT, on the other hand, the source flux can be either specified or it can be read in by the program, from a file produced by

CIGARET.

It is equally necessary to specify the heat losses. If the existence of the cigarette is neglected for the

moment, so that the view factors need not be considered, then

<t>.«r = ''(^,-3;)-€or/ (4)

where the first term is the convective cooling term and the second is the radiative cooling term. In

contrast to the source term, these losses are explicitly included in the program. Here, the dependence

of Tg on X, y, and t has not been shown explicitly, e is the emissivity of the surface and a is the StefanBoltzmann constant. (Note: we have used e where sometimes a, the absorptivity of the surface, should

be used. When there is thermal equilibrium between the radiation field and the hot surface, then a =

e. We have assumed the latter, for simplicity). Evidently, the loss rate from the cushion surface

increases as it heats up. h is determined by the laws of fluid flow, as well as the thermal properties of

air; the simplest way to find h, however, is to use well-known expressions, usually derived from

correlations. This will be further discussed below.

Along with the partial difl'erential equation (12), we must have boundary conditions: for the top surface,

the net flux entering the surface at each point is connected to the net flux according to

*,.(-.y.o = -Kf^^^^] z = o. t>o (5)

where <^>net = <^in - <^out • (6)

When the cigarette lies on the substrate, the mean convective heating flux over the heating region is given

by

If hq is the heat transfer coefllicient for quiescent air, then the surface heat loss for the areas away from

the cigarette is
(8)

The boundary condition is greatly simplified if we can say that equation (8) is valid over the entire

surface. That is readily achieved by using the model convective heating flux <f>^*,

*; = ^.(Tli, - T,) - h^iT, - TJ (9)

hjn is found as follows: the mean initial convective flux from a cigarette was measured to be 37 kW/m^.

TaJcing the mean value of the cigarette surface temperature to be about <T>cig « 450 °C, the mean

heat transfer coeflicient in this area is hj^ « (37000/430) = 86 W/m^ K. h^ is found from Table 7-1

of Holman (1981): the Nusselt number for a heated, horizontal, upward-facing surface is

where Ra is the Rayleigh number,

surface, we find:

Nu = 0.54 (/2fl)*/* (10)

Since h = kNu/L, where L is the characteristic dimension of the

h = 0.54 k

1/4 (11) a vL

where /S is the volumetric coefllicient of (gas) expansion. Taking the value for T^ = 500 K « 227 °C

as an approximate mean temperature of the surface, and L = 7 mm, we find:

hq = 9.7 W/m^-K

For the other surfaces, the boundary conditions (b.c.) can be expressed in diff'erent forms, depending on

whether the slab for which we are making calculations is exposed to the air, or is embedded within the

cushion. If the slab were the entire cushion - that is, the other five surfaces are exposed to the air, then

the b.c. would be that the heat fluxes KdT/dx, KdT/dy, KdT/dz, are given by terms of exactly the same

form as (4), except that h is difl'erent for the vertical and the downward-facing horizontal surface. On

the other hand, the numerical calculation we are using makes using the entire cushion prohibitively large.

We therefore consider a subsection abstracted from the whole cushion, and the other five faces are surrounded by more foam. If that foam were a perfectly insulating material, then the b.c. at those faces

is the adiabatic condition; that is, no flux crosses the surfaces: KdT/dx = KdT/dy = KdT/dz = 0. This

eliminates any heat losses from the slab, and makes the calculated temperatures somewhat higher than

they really are. A simpler b.c. is to assume that the temperature at the five faces is constant and

unchanging, remaining at the original (usually ambient) temperature. This produces a larger temperature

gradient than actually prevails, and errs in the opposite direction. That is the default b.c. used in the

calculation, but which of the two b.c.'s is used is determined by an input parameter.

2. Heat Transfer

The partial difl'erential equation which describes the conduction of heat in a solid, when there is no

radiation heat transfer, does so by giving the rate of change of temperature at every point within the solid.

It is (see, e.g., Carslaw and Jaeger (1959))

pM! = div(K gradr) + 5 (12)

6
where T, which is a function of position and of time, is the variable for which we are trying to find a

solution. The other symbols in this equation are: p is the density of the solid, c its specific heat, k its thermal conductivity, and S is any volumetric heat source or sink, p, c, k, and S may vary with position

and with T. Because these may vary in equation (12), it is conceptually simple to consider a substrate

which consists of layers of material. Indead, it will be straightforward to implement this in practice, as well. Note that because we only have derivatives of T, T need not be the absolute temperature: it may

be taken to be that relative to a convenient reference temperature. In other places, such as in equation

(2), it must be taken to be absolute {i.e., in Kelvins).

The first term on the right-hand side describes the diffusion of heat in the solid. If there is any release

of heat ~ usually by combustion - at the (interior) point in question, it is given by S(x,y,z,t). A general

expression for S is

5 = E^^c.^^v-|^ (13)

where R; is the reaction rate (in kg/m^-s) and H^,j is the heat of combustion (in J/kg), for the ith reaction. If there is endothermic pyrolysis or evaporation, we have the last term on the right-hand side

as well, where is the latent heat of evaporation or pyrolysis. An explicit expression for Rj will be

given in Section II. B. 4.

Since the furniture (apart from the frame) consists of a fabric-covered pad, it is clear that the program

must take at least two layers (with diff'erent properties) into account. Therefore the program was written so as to permit diff'erent values for the relevant thermophysical constants p, c, and k in each layer. In

fact, generally there is not perfectly intimate thermal contact between the fabric covering and the padding;

there is a small but sometimes significant intervening air gap. Normally, one would place a node within

this gap, in order to take a third layer into account; because of the thinness of the gap, and other technical

difficulties, however, a diff'erent treatment of the eff'ect of this air gap has been devised: the gap can be

represented in terms of its "thermal resistance." (See Section II. C. 3.) This has been programmed and

successfully tested.

In writing equation (12), we have made the simplifying assumption that the cushion is totally opaque; that

is, there is no radiative transfer of heat through the cushion. If there were, the equation describing heat

transfer through the solid would become still more complicated. However, the fabric and foam can each

be thought of as porous, consisting of solid parts interspersed with void spaces. Then taking forward

radiation transfer in those spaces, it is possible to incorporate a first approximation to (one-d) radiation

transfer, as shown by Kunii (1961)

where $ is the void fraction, and and are the solid and gas-phase thermal conductivities,

respectively; they are each functions of T. Dp is the mean pore diameter, and

h, = 4€aT^ (15)

Equation (14) was equation (5-50) in Gann et al. (1988). Finally, we must note that at a giv^

temperature, the thermal conduaivity is proportional to the density (also see Section n.C.5):

*(P) = (p/Pa)*(Pa) W

7
3. Time to Ignition

An experimental test for the time to ignition is described in Section II.D. We can readily make an orderof-magnitude estimate for this time, based on equation (12). For the very special case that the substrate

is initially at the uniform temperature T^, that it is homogeneous, isotropic, inert, and semi-infinite, that the thermophysical characteristics p, c, and k of the material are independent of the temperature, and that the problem is one-dimensional (i.e., the incident flux is the same everywhere on the plane z = 0,

resulting in the slab being subjected to the uniform net heating flux <^>net(t))» then equation (12) can be

solved analytically and explicitly, and it can be shown (Carslaw and Jaeger, 1959, p. 76) that the surface

temperature is given by

r.(o = r..-J^/*=^' (.7)

For the still simpler (unrealistic) case <|>J^^ = constant, this can be immediately integrated, and we obtain

TX,)-T*^*^ (18)

When the left-hand side reaches the ignition temperature Tjg, the time elapsed must be the time in the

right-hand side of the equation, and may be called the ignition time:

net )

(19)

Although not a single one of the simplifying assumptions required to obtain equation (19) is valid, this expression is nevertheless very useful as a guide to the form of the dependence of tj^. In particular, we

will use equation (19) in order to make sense out of the experimental observations, m Section n.E.l.

[We can, at some cost in computation time, get a somewhat more realistic result by relaxing one of the

assumptions above: it is possible to write, starting with equation (17), an explicit expression for T5(t) for

the case where the net flux is not constant, but results from some impressed external flux, diminished

by a Newtonian (/.e., linear) cooling loss; that is, where we can write ^^^^ = h*(Tg - Tg); see Quintiere

(1988)].

4. Pyrolysis of the Substrate

As described above, we assume that the substrate consists of a thin fabric covering a relatively thick foam

pad, with a very thin air gap between them. Each of the two materials will be heated and can pyrolyze.

In this Section we describe how these reactions are calculated. Pyrolytic reactions can generally be

expressed in the Arrhenius form

R,iZ P,. P,) = A p7 p:, e -^^"^ kg/m' s (20)

8
where Pf is the density of the fuel and oxygen; E^^^ is the activation energy for the reaction,

and R is the universal gas constant (not to be confused with tiie reaction rate Rp in equation (20)). It is sometimes convenient to define an "activation temperature" T^:

= EJR (21)

In the remainder of the paper, R is used only to represent a reaction rate. All the quantities in equation

(20) except the pre-exponential factor A are at the location (x,y,z) within the porous solid and at the time

t; this dependence is not shown explicitly in order to minimize the complexity of the expression.

Similarly, the subscript i, corresponding to ith "reaction" has been suppressed from R, A, pf, m, n, and

E^. If the reactions were truly elemental Arrhenius reactions, then m = n = 1. However, equation (20)

is really a model equation, and therefore m and n are purely empirical parameters (just as A and E^^ are),

and need not be integers.

Oxidative as well as non-oxidative pyrolysis (thermal decomposition) of materials have been included..

The former reactions are generally exothermic, while the latter are generally endothermic. There may

be several oxidative reactions as well as several non-oxidative ones; char oxidation, when it takes place,

is the last step. In this model a maximum of three reactions is permitted for each material. The oxidative

reactions are exothermic; we assume that the rate of each step (reaction) is adequately described by an

Arrhenius equation of the form (20).

It had been hoped that the diflFerent pyrolysis steps take place at sufficiently different temperatures, that the reaction sequence for the ignition of the fabric could be taken to consist of a small number of steps

which follow each other sequentially. That would correspond to the equations taking the form

^ = E = E (P - P. Pox ^^(-TJ T) (22)

i i

where pj is the density of the fuel at the end of the i** reaction step. When the reaction produces a gaseous species which escapes, the reaction rate is given by

/? = --^(23) ^ dt

More generally, if species i is reacting and producing species j, then

^ = -i? = (24)

dt ' dt

We will here discuss only the pyrolysis of the fabric, which is of principal interest: when smoldering

ignition is produced, it is almost always first initiated in the fabric. Pyrolysis of the foam is discussed

briefly in Section n.D.2. The fabric of particular interest is cotton duck; cotton is principally cellulose.

Analysis of the reaction kinetics of a cellulosic paper by Kashiwagi and Nambu (1992) showed that it could be described by three reactions, all of which will proceed simultaneously (though obviously at diff'erent rates). That means that the equations describing the creation and destruction of different fuel species must be included. Virgin material goes to char, via two pathways: degradation and oxidative

pyrolysis. Each gram produces n^. grams of char. The char then goes to ash, via char oxidation. Each

gram of char produces n^ grams of ash. Thus one gram of material produces n^n^. grams of ash. The

governing equations are

9
and

The density of the solids is

f>y = -Rd-Rop (25a)

Pa = ^aK (25c)

Pj = Pv Pe * Pa (25d)

Hence

p, = -{l-n^){R,-R^)-{\-n^)R,, (25e)

The subscripts denote:

a = ash op = oxidative pyrolysis

c = char s = solid

CO = char oxidation v = virgin material

d = degradation

Besides the surface heat source described in Section II.B.l, there are also internal heat sources (and

sinks), given by equation (13). The heats of combustion must be supplied by the user of the model; those

for cotton are given in Section II.D. The global kinetic constants given by Kashiwagi and Nambu are

also listed in Section II.D. If the reaction rate of one or more exothermic reactions becomes high, for

any reason, the rate at which energy is being generated will exceed the rate at which it is lost, and a "thermal runaway" will ensue, very similar to a chain reaction. During the runaway, the rate of pyrolysis

(and of heating) is limited by the availability of oxygen. Since O2 has 0.001 times the density of the local

fuel, it would immediately be exhausted, but for the supply which (a) diffuses in from adjoining cells,

and that which (b) diffuses in (or is convected in) from the surface.

5. Gas Diffusion

As oxygen is depleted in the regions where combustion (oxidative pyrolysis, char oxidation) is taking

place, the concentration gradient which results will induce diffusion of oxygen into those regions, from

the surface as well as from adjacent, oxygen-rich regions, assuming the medium is porous and permits

diffusion. Similarly, the gaseous (and other) products which are generated - mainly CO2 -- build up in local concentration, and diffuse away.

The equations which describe the rate of change of species concentration are analogous to equation (12)

for the diffusion of heat; assuming no convection, the ith species density is governed by

= div(Z)^gradp,)+5, (26)

where Dj^ is the diffusion coefficient for species i in the background "o," and Sj is the source/sink term.

If the diffrision coefficient increases with T, as it normally does, then the CO2 diffusion, for example,

will take place preferentially towards the hot regions, i.e. , towards the surface. The migration of gases

is expected to have a minor effect on the heat transfer, and hence on the temperature distribution, from

10
this effect. However, the diffusion of oxygen is important: the rate at which oxygen enters the reacting

region sets the limit on the reaction rate. A porous medium is one consisting of solid particles embedded in a gaseous medium. Put in a different

way, a gas molecule cannot penetrate any of the solid phase, but is able to traverse the entire medium,

either because of a pressure difference, or merely from the "random walk" which constitutes diffusion.

The diffusivity of a gas species in such a medium is really the diffiisivity of those gas molecules through

a gaseous medium consisting of the same or some other gas, where the solid volumes are excluded. That

is the effective diffusivity.

It has been found experimentally (Szekely et al., 1976) that, approximately,

D,ff = Z),*" (27)

where $ is the fractional void space; this is also referred to as the "porosity" in the literature. The

relationship is strongly dependent on the structure of the (granular) material; thus. Figure 2.4 from

Szekely et al., 1976, shows that for particles of mica,

DeffID, = (27a)

whereas for sand, a bed of glass spheres, carborundum powder, and table salt,

^effl^c = 0.677*''* (27b)

for * < 0.7. As shown in Appendix E, $ « 0.6 for the cotton fabric that will be our main focus of

interest; hence equation (27b) is the relationship to be used.

Here is the diffusivity of O2 in N2. From equation (16.3-1) of Bird et al. (1960), we find that for

O2 in N2, the temperature dependence is

/ 'p \1.823

D (J) = 0.199 —-— cm^/s (28)

1,293.16 j

6. Boundary Conditions for Gases

There are two cases of interest: first, when the air above the substrate is quiescent, and second, when ~

as was the case in the experiment to be described in Section II. D -- there is a stream of air impinging on

the substrate. We first obtain a general expression:

When the air above the substrate is quiescent, and there are oxidative reactions taking place in the

substrate, the concentration of oxygen molecules in the air decreases towards the surface, through a boundary layer. The maximum possible reaction rate is obtained by assuming y^ = 0, as can be seen from equation (29). On the other hand, this is clearly not possible, since we must have a finite concentration at the surface in order to get any reactions at all. It is possible to obtain a value for yg analytically, if one makes the simplifying assumptions of an isothermal substrate and a first-order reaction

rate (Ohlemiller, 1991).

11
Mass transfer is a (molecular) transport phenomenon, just as is heat (momentum) transfer. Thus,

analogous to the first term on the right-hand side of equation (4), the mass transfer of oxygen across the

surface of the solid can be written as

(29)

where is the ambient oxygen mass fraction and Yg that at the solid surface; the x,y dependence of the

latter has been put in explicitly in order to emphasize the origin of the spatial dependence of iho2"- This

is the boundary condition.

Note: The velocity at which oxygen enters the surface is

V, = =y(Ya-Yi) (30)

(see Gann et al., 1988, pp. 159-160) where D is the diffusion coefficient inside the solid. Thus

(31)

where y is the normalized mass transfer coefficient for oxygen traversing the boundary layer above the

substrate. In SUBSTRAT (see Gann et al. (1988), equation (5-64) (from Muramatsu, 1981)), the

following equation was used for y:

Yj, = 6MxlO

7-2.75(7, - rj(7 + 123.6)

V4

cm/s (32)

where the subscript "b" refers to "boundary layer". On the other hand, in this paper we obtain 7 in a

different way; we will afterwards use equation (32) in order to compare the results. The way 7 is obtained is as follows: Because of the similar origin of mass and heat transfer, one can often use the

Reynolds-Colbum analogy. It is not difficult to show, from the treatment in Chapter 3 of Treybal (1955),

that the Reynolds-Colbum analogy leads to

Y

= m/s (33)

where h is the heat transfer coefficient appropriate to the problem, Pr is the Prandtl number, Sc is the

Schmidt number (Sc = vfD), and y is the kinematic viscosity. Once we have 7,,, equation (31) gives

kgas, to be used in equation (29).

Case A. Ouiescent air

Here we have

h = kNuH^ (34)

where k is the thermal conductivity, 1^ is a characteristic dimension for the problem, and Nu is the

Nusselt number. We can thus rewrite equation (33) as

12
(35)

In quiescent air, there is a constant movement of molecules in all directions; the flow in any one direction

is exactly compensated by the flow in the opposite direction, normally. If, however, the gas is near a boundary which "absorbs" some of these molecules (as is the case for oxygen impinging on a reactive

substrate) then evidently the return flow is smaller, and hence the impinging flow is not completely

nullified. The result is that there is a gradual decrease in the concentration of those molecules as the

surface is approached. This region is referred to as the "boundary layer."

The quantities in equation (35) are to be evaluated at some characteristic (mean) temperature of the

boundary layer. That temperature can be taken to be the mean between T^ and Tg. For quiescent air,

with the surface horizontal and facing upward, is the boundary layer thickness. The expression for the

Nusselt number is given in Section n.B.l.

Case B. Impinging Air

The Nusselt number for this case is given by equation (DIO), Appendix D. Inserting that into equation

(35), we obtain

where a is the thermal diff\isivity of air. Re = wjjv is the Reynolds number and u^ is the incoming

velocity of the impinging flow. We now show that when oxygen consumption is significant, it takes place in a very narrow layer; this

is analogous to the flame sheet approximation for combustion in air. As a result, we may take the

concentration of oxygen in that layer to be some appropriate average value, rather than having to actually

program in the species diffiision equations, via equation (26).

In this Section, it is shown that a reasonable approximation to the overall reaction rate can be obtained

without explicitly solving the gas diffiision equation, equation (26), for the difl\ision of oxygen in the

substrate. This is the approximation made in the model.

In the experiment to be discussed in Section II. D, a jet of air is directed downward at the substrate. Hence there is no boundary layer in this case, and the oxygen concentration just outside the surface is 21%. When the reaction rate is low, the oxygen concentration in the top of the substrate will not be

much aff'ected by the slow reactions. When the reaction rate becomes high , on the other hand, the mean

oxygen concentration in that thin layer becomes low. It cannot get too low, however, since (as discussed

in Section II. B. 6), the reaction rate would then fall. Hence a quasi-steady state is established. Thus,

in the "runaway" phase, the reaction rate must be given by the boundary condition equation (29).

We will now argue that when significant oxidative pyrolysis is taking place,

- 0.767 -v/^Pr^ *" 5c -2/3 (36)

7. Oxygen Diffusion

13
(a) taking a constant value for the oxygen concentration ([O2]) in the substrate is a valid

approximation.

(b) most of the pyrolysis occurs in the surface layer, and therefore

(c) we obtain an adequate approximation to the pyrolysis rate, without explicitly including

the oxygen diffusion equations.

When the surface temperature is relatively low, there is negligible pyrolysis, and the oxygen distribution

[02](x) in the substrate is approximately uniform. When has risen to the point that perceptible

pyrolysis is taking place at and near the surface, the [O2] profile in the substrate dips as the surface is approached from below (i.e., from within);it must reach a minimum and then rise again at the surface,

because of the diffusion of oxygen from outside. When the temperature reaches values so high that [O2]

is pulled down to negligible values near (but not at) the surface, the principal source of oxygen in the

reaction zone is from the air diffusing in from the surface; the region where there is significant oxidative

reaction is then highly localized near the surface, with [O2] being highest at the surface, and falling

rapidly until it reaches negligible values. The characteristic oxygen penetration distance is 6, and it is clear that 6 is a steeply falling function of T. [At some further distance in, the oxygen must diffuse

towards the surface from other (deeper, or peripheral) parts of the porous substrate, because of the

concentration gradient, so that the profile must rise again as we go deeper. It can do so because at that depth the temperature has fallen sufficiently to "freeze out" the (oxidative) reaction rate, i.e., it is negligibly small].

Suppose that is high enough to drive [O2] to near-zero at some depth. Then we may write, crudely,

Yix) - Y^e-''"' (37)

with Yg = surface value of [O2], and 5 = characteristic penetration depth. If we write the reaction rate

in the form

R[x,T(x)] = Y(xyF[x,Tixy\ = Y^Gix) (38)

then the total reaction rate per unit area is

to" = rY(xTG{x)dx (39)

J o where G(x) is a fairly steep (decreasing) function of x. We may thus further write (still very crudely)

Then

/ 1 \

J_ JL

1/2 m" ~ (41) fy:e'"'l''G,e-^''^dx = G^^

We use the mean-value theorem, and write the integral in equation (41) as

("rG^e-'^'^dx = Y'G^Q^n (42)

J o Then

14
nQ-'X"^ (43)

1 + —

We now need an estimate for Bib. Clearly, the reaction rate must be significant for a distance comparable

to the penetration depth 6: else, if the reactions "froze" before [OJ becomes negligible, [O2] would not

fall below some appreciable value. Moreover, the reaction rate cannot be significant below that depth,

precisely because there is so little oxygen below. Thus ^ = 6. For n, we take n^, = 0.5, as given in

Kashiwagi and Nambu (1992). Then equation (43) yields

y = 2y^/3 (44)

Next, consider the efl'ect of discretizing the equations (Section II. C). If 6 > Ax, then assuming that the

reaction takes place in the top cells only, clearly underestimates the reaction rate. However, when 6 is large, the total reaction rate is very low. That is, while the fractional error is large, the absolute error

is small, and (we shall show) typically negligible.

Assume, for the sake of simplicity, only one oxidation reaction. Assume, also, that the temperature and

reaction rate are uniform in a thin sheet of depth 6 at the surface, where T is maximum. (This is essentially what is implicitly assumed in the numerical calculation.) Then if the stoichiometric

fuel/oxygen ratio is r, the fuel bumup rate is Rb =

m'l = rm;;^ (45)

where R is the reaction rate (g/cm^s or kg/m^s). If 6 is greater than the thickness of the surface cell, then

evidently the reaction rate in the top layer is limited to RAz; whereas if 6 < Az, equation (45) gives the

limit. Thus the reaction rate in the top cell is

RjopAz = min[RAz, rm^,^] g/cm^-s (46)

Note that when the reaction rate is so high that 5 < Az, the concentration at the bottom face of the top

cell(s) is close to zero; therefore a reasonable value to take for [O2] in the cell is on the order of the mean

value, 11%, (or 0.12 for the mass fraction); or, according to equation (44), <Y> « %(0.2318) = 0.15.

We now make a calculation to find a typical value for 6. As is clear from equation (45), we must begin

by calculating a reaction rate R. In order to find R, we must anticipate some results from Section II.D.2.

Using equations (79) and (80), we find that

Rop = Pvkop = Pv l-SxlO^^ Y<^-5 (WdAVJ»-3 exp(-160/RT)

and

Rco = Pckco = Pc 3.4x10^1 Y^'-^* (W,/W„) exp(-160/RT)

The meaning of the subscripts is made clear in the list below equation (25e). These rates are in kg/m^

min. We must choose characteristic values for Py, Y, W^jAVq, p^,, W^/W^, and T, in order to get

estimates. Reasonable values are: Wj = Wo/2, W^ = Wo/4, and Y « 0.11

For Py and p^, we assume that the volumes are unchanged, so that

15
Pv/Po = WdAV, and p^o = W,AV,

Finally, we take = 1560 kg/m^. We then find that

k^p = 2.02x10^^ exp(-160/RT) min'*

and

k^o = 1.52x10*° exp(-160/RT) min*

The universal gas constant R = 8.31447 J/mol-K, so that the activation temperature for both reactions

is

Ta = 160,000/R = 19244 K.

For a temperature in the vicinity of the ignition temperature, say, 400 °C = 673.16 K, we then have

Kp = "^-^^ min* and k^^ = 0.00584 min"*,

yielding

Rop = PXp = 101 itg/m^ s

and Rco = Pc^^co = 0 038 kg/m^ s.

Thus at 400 °C, with these assumptions for and p^, oxidative pyrolysis takes place more than 2600

times faster than does char oxidation, and we may therefore take

R = Rop = 101 kg/m^ s

Next, we need to have the rate at which oxygen enters the medium, m^^ . Equations (29) and (31)

give the required value. For the ignition experiment, the appropriate expression to use to find y is equation (36). The temperature at which the various quantities must be evaluated is that of the purging

gas. Examination of Table D-1 in Appendix D shows that Tg is quite high. Assume T- = 800 K.

Equation (28) yields 0^(800) = 1.24 cm^/s. Then the Schmidt number is Sc = 0.663, and

7 « 8.10x10"^ V^e li^ m/s

(with in meters). The characteristic air velocity for calculating Re is given by eauation (D15) in

Appendix D. Finally, the characteristic length is the standoff distance, 6. With 6 = 5.4 mm, we find

Re = 11.06 and y = 4.99 cm/s. (Incidentally, if we had quiescent air at 800 K, Muaramatsu's

expression, equation (32), yields 7^ « 2.35 cm/s).

Finally, taking = 0.232 and = 0, equations (29) and (31) then yield

mj; « 13.6g/m2s

We must also have r (the stoichiometric fuel/oxygen ratio) in order to use equation (45) to get the

penetration depth 5. The stoichiometric air/fuel mass ratio for wood is S = 5.78, close enough to that for cellulose. Hence r « 0.746, and we find 6 = 0. 10 mm. Thus the penetration depth 6 is indeed

much smaller than the layer thickness Ax (which is of the order 0.5 mm).

16
For = 300 °C, on the other hand, we find k^p = 0.053 min"', so that = 1.38 kg/m^ s. Then

6 « 7.3 mm, several times Ax, but the resulting mass-loss rate is only 0.04i g/m^ s, which is indeed

negligible.

For the case of quiescent air above the substrate, the mean value for the oxygen concentration to be used

in the top cell during the runaway, is midway between the value at the top surface (i.e., at the bottom

of the boundary layer), and at the bottom of the top cell (presumably, close to zero).

The only time that we cannot justifiably make the simplification that the oxygen concentration in the top

layer is either that at the top surface or half (or 2/3) that value (during runaway), is when the reaction

rate is intermediate between the very low values and the runaway value. This period should be relatively

short, and the error introduced by tiiese simplifications should not be large.

C. NUMERICS

In general, it is not possible to solve equation (12) analytically, in spite of the simplifying assumption of no radiation heat transfer, and even for the case S=0. That is, to write down an explicit expression

which gives T in terms of the inputs. This is so, because of the nonlinear and nonuniform boundary

conditions. It is therefore necessary to resort to a numerical procedure, which is that incorporated in TMPSUB2 and SUBSTRAT. Hereafter, we shall only refer to the latter, when something applies to both,

as is the case here).

1. Introduction

The development of TMPSUB2 centers around the capability to simulate transient

heat transfer. A one-dimensional heat conduction problem provides the simplest example to illustrate transient simulation methods.

Figure 1 shows a portion of a one-dimensional conduction problem in which the

material has been divided into thin layers.

This example will be described by physical

instead of mathematical arguments following

the description given by Clausing (1969,

pp. 157-213).

The figure focuses on a representative layer

of material of thickness Ax; centered in a node at coordinate Xj. This material layer is represented by a single temperature Tj and a corresponding thermal conductivity density p ^, and specific heat C;. Assume this layer has a surface

area of magnitude A in the Y-Z plane. The distance between nodes i and i-1 is given by Sxj (= Xj - \{).

Subscripts i and n refer to positions in space and time, respectively. Subscript n is temporarily

suppressed until it becomes necessary to consider time, in the following equations.

The instantaneous internal energy of layer i is given by

= pcAVT. = pc^x.AT, = C^T^ (47)

layer I -AX,- —•

-• X

: ^

1 X1+1

Figure 1. One-dimensional conduction

17
where T is the absolute temperature. This equation also serves to define the heat capacity, C, assigned

to node i. Heat is transferred to and from layer i by three methods:

(1) conduction from layer i-1 at the rate

qi. = {T,_,-T^)K^_Albx, = K,_iT^.,-T.) (48)

(2) conduction from layer i + 1 at the rate

= - 7',) K,,^ / 6^,,i = K,^{T,,, - T,) (49)

(3) internal heat generation or radiation absorbed within the layer, S j.

Equations (48) and (49) define the thermal conductance, K, in each direction. The possibility of a thermal conductivity which varies with position is explicitly taken into account here using labels + and —: K^__ and /cj^.. A good first order approximation for K^_ is the harmonic average of /<j and /Cj_i:

K,. . (50)

with a similar expression for

The change in internal energy of this layer between time tn and time t^+j is given to first order in At by

AU, , sC ,r. Z = At(Q. +Q. +q.) (51) i,n-n+\ i,B+l i,n i,H "•^V"*,. i,+ i,/ v /

where the time dependence has now been put in explicitly. There are several common solutions to equation (51) depending on when the heat gains are evaluated. One solution involves evaluating at time n

c,„.iT,„^i = c,„r,„ . Ar[^,.(7;..,„-7;.„)./:,,(7;,j„-7;.„).5,j (52)

which is the standard Euler explicit time integration formula. "Explicit" means that T; „+i can be

directly computed from values known at time n. On the other hand, evaluating at time n+ 1 gives

which is Euler's standard implicit time integration formula. "Implicit" means that Tj jj ^. j is computed

from other values also evaluated at time n+ 1. These values depend implicitly on each other and must

be computed by a solution of simultaneous equations.

Clausing (1969, p. 190) also gives a discussion of stability in terms of thermodynamic laws. Rearranging

equation (52) to solve for Tj gives

where the time subscript, n, has been added to the K terms to indicate exactly when these values are evaluated. For the sake of argument, assume Cj = Cj „ (=0; ) for simplicity. There is no solution

if Cj = 0. If Cj is sufficiently small or At sufficiently large, then ( Ki_ „ + Kj^ „ ) At / C; > 1, and

18
as Tj n increases Tj j^+j must decrease, and vice versa. This is thermodynamically impossible. It shows

up in a numerical solution as oscillations, i.e., "instability," in the node temperatures at each time step.

These oscillations tend to quickly increase to totally meaningless values. In general, the smaller the

thermal mass of the element, the smaller the time step needed for a stable explicit solution. This suggests

a simple technique to determine the minimum stable time step for any element in the system.

Thus, rearranging the implicit (53) to solve for Tj n+i gives

J. ^

^i.n^i.n ^ ^ (^i-,n*l '^i-l.n+l '*'^i*l.n*l ^i*l.n*l ^^i.n+O ^55-)

This equation shows none of the computational or thermodynamic problems of equation (54), indicating

that the standard implicit method is stable for all time steps.

The spatial discretization error (for a uniform grid) for the standard explicit and standard implicit methods

is proportional to (Ax)^ (for a variable grid, the accuracy is reduced somewhat; see Section II.C.4). The

time discretization error is proportional to At.

The standard explicit and standard implicit methods err in opposite directions. Therefore, a more

accurate solution can be obtained by combining the two methods. Expressing this combination generally

in terms of a parameter /3 gives

where 0^/3^1.

/S = 0 corresponds to the standard explicit method,

/S = 1/2 corresponds to the Crank-Nicholson method,

/S = 2/3 corresponds to the Galerkin method, and

/3 = 1 corresponds to the standard implicit method.

For jS ^ 1/2 this method is unconditionally stable, although the solution may be oscillatory. For /3 >

3/4 (approximately) the solution is stable and non-oscillatory. For /S = 1/2, the time discretization error

is proportional to (At)^.

The methods presented above extend directly into three dimensions. For a Cartesian coordinate system

and a cell of dimensions Ax; by Ayj by Az^, equation (51) can be rewritten to account for conduction

from the six adjacent cells in the 3-D system:

where

19
= ( '^ij,k-x - Tij,t) ^k- ^yj i^h

<ik. = ( '^ijMx - '^ij.i) I *^)t.i

and Sj

j k represents other heat added directly to the cell.

2. Boundary Conditions

Special treatment is required for cells on the boundaries of the region being modeled. In particular, the

nodes which represent the cells are placed on the boundary, rather than at the center. Referring to the

one-dimensionaJ example in Figure 1, the surface layer is only half as thick as the others.

An adiabatic boundary condition (b.c.) is handled by setting the appropriate heat flux terms in equation

(45) to zero. A constant temperature (isothermal) b.c. is handled by leaving the temperature unchanged.

The surface of the substrate (z=0 plane; see Figure 3) transfers heat to the environment by convection

and by radiation. The convective heat gain for cell i,j,l is given by

(5c)i.y.i =(';-7',;.i)AAx,Ay, (58)

where T^ is the temperature of the surrounding (ambient) air, and h is the heat transfer coefficient,

A positive value of S^, represents a heat gain by the cell. The radiative heat gain is given by

where

Ta is the temperature of the surrounding surfaces, a is the Stefan-Boltzmann constant, and

£ is the emissivity of the fabric (assuming that a, the absorptivity, = c).

The temperatures in equation (59) must be absolute (Kelvin) temperatures. In TMPSUB2 the temperature

of the surrounding surfaces is assumed equal to the air temperature. (Note that equation (59) can be

rewritten in an apparently linear form similar to equation (58):

This form will be useful further on.

The heat flux from a smoldering cigarette to the fabric is, for the "prescribed flux" choice of input,

represented by the following pair of equations:

{x^x+vt) (60a)

20
ix<x^+vt) (60b)

This heat flux is converted to a heat gain, Sg, by multiplying the flux at the position of the node by the

cell surface area. These values are adjusted slightly so that their (discrete) sum is equal to the total heat

flux represented by equations (60a,b):

3. Air Gap

Since the furniture (apart from the frame) consists of fabric-covered padding, it is clear that the program

must take at least two layers (with diff'erent properties) into account. Therefore the program was written so as to permit diff'erent values for the relevant thermophysical constants p, c, and k at each node. In

fact, generally there is not perfectly intimate thermal contact between the fabric covering and the padding:

there is a small but sometimes significant intervening air gap. Normally, one would place a node within

this gap, in order to take a third layer into account; because of the thinness of the gap, and other technical

difficulties, however, a different treatment of the efl'ect of this air gap has been devised: the gap can be

represented in terms of its "thermal resistance."

The optional air gap between the fabric and

the padding is modeled by assuming one- dimensional heat transfer. Figure 2 shows

the basic configuration and nomenclature for an air gap of thickness s between cells i,j,k

and i,j,k-t- 1. s is small relative to the separation between the cells, Although

pyrolysis of the fabric and the possible

melting and/or pyrolysis of the padding may

well change the dimensions of the air gap,

the simplifying assumption is nevertheless

made here that the air gap has fixed and

uniform dimensions. The heat transfer

between these cells is given by

(61)

where K is an implicit function of qij,k+' ^ij.ic ^i,j,k:+-

The overall heat transfer coefficient is given by

and

<|),(x,)'.f) = 4>;,«exp

MuiBOill

IPH

^.

K,( ) Ko

padding C9

Figure 2. Air gap heat transfer model

21
K =

1 1 1 + — +

where the fabric conductance is

the padding conductance is = lLx^Ly.yL^j^._^^l S^^j ,

the fabric (bottom) surface temperature is ~ ^ijk~ ^ij t+ / ^ij k '

the padding surface temperature is = T^j^^^ + q^j / k.^. ^.^j ,

the radiant conductance is

and the convective conductance is

The conductance of the air, /tajr, is evaluated at the average of Tf and Tp. These equations are solved to give K and qij,k+ during the overall process used to compute cell temperatures.

4. Variable Grid

The heat from the cigarette spreads into the substrate by conduction. In order to correctly estimate the

temperatures near the peak, it is important that conduction to the outer boundaries (at x=0, x=x^ax»

y~yinax' ^d z=z^ax) ^® negligible, since it cannot be known a priori. This requires a relatively large

region relative to the size of the cigarette heat flux pattern. However, setting Ax, Ay, and Az suflliciently

small to achieve the desired accuracy in the conduction calculation can result in an extremely large

number of cells ( N = (x^^g^/Ax) • (y^^/Ay) • (z^^/Az) ) with correspondingly large memory requirements and computation time.

The heating flux from a smoldering cigarette rises from negligible values to a high peak, on the order

of 60 kW/m^ over a region only a few millimeters in extent. In order to follow this faithfully, the region

must be covered by a mesh which is fine enough so that there are no changes from one mesh point to another large enough to produce numerical inaccuracies or instabilities. Thus, the required size of the

grid is inversely proportional to the temperature gradient, where steep temperature gradients occur only

near the point of peak heat flux. Therefore, a variable grid is used. This grid consists of a few constant- width cells near the peak followed by cells of regularly increasing size to the outer boundaries. The

increase in cell size is based on geometric progression, and can be diff'erent for each axis. Thus, for

example, Ayj+i = Ry Ayj (R^ ^ 1). The general efl'ect of this variable grid is shown in Figure 3. This

is easy to implement m equation (45) which explicitly incorporates the grid sizes. The variable grid gives

22
results which are less accurate than a constant

grid. A benchmark test indicates that the results for R = 1.23 differ by < 0.11% at t = 100 s. from those obtained for the constant grid case, with still smaller errors for R closer to unity.

(See Table A-1, Appendix A.)

Note: In choosing a grid, the user determines four

items, in each coordinate direction: the width of

the constant-width cells in the fine-grid region,

the number of such cells, the total number of

nodes along that axis, and the total width in that direction. The program then does the arithmetic,

and finds the corresponding value of Rj. If that is < 1, an error message will be returned. Likewise, if too many grid points result from the

attempted selection, an error message will result.

The point of peak heat flux moves as the cigarette smolders. TMPSUB2 adjusts the x-coordinates to keep

the fine grid region centered on that peak. This adjustment is made by

(1) computing the new x-coordinates of the shifted grid,

(2) computing a cubic spline curve for the temperature in each row of cells in the old grid,

and

(3) using the curve to compute the cell temperatures in each row at the new grid positions.

The standard explicit algorithm for solving equation (12), without pyrolysis included, was checked in

several ways, principally by comparing its predictions against known analytical solutions. (See Appendix

A.) These checks showed that the numerical procedure, and the computer program for implementing it, are correct, effective, and accurate.

5. Variable Thermal Properties

The conductivity and specific heat of the substrate materials are known to vary with temperature.

TMPSUB2 allows the user to describe this variation. The values entered for k and c are converted to cubic spline curves giving x(T) and c(T). The conductivity is further adjusted within the program to take

into account the fact that the thermal conductivity is proportional to the density:

/c = p k(T) / (62)

where is the original density (recall that the assumption of constant k, p, and c, in Section II. B. 3, was

a special case used only in order to obtain equations (17) to (19)).

6. Pyrolysis

Consider pyrolysis from the finite difference viewpoint. We begin with a cell having constant volume,

AV (= Ax Ay Az) and temporarily ignore the subscripts i,j,k. This cell contains a mass Py AV of virgin

(unpyrolyzed) material, a mass p^. AV of char, and a mass p^ AV of ash. The numerical method must

V Z

Figure 3. Substrate coordinate system

23
keep track of the density of each material as a function of time; in this program, this is done via the

difference equation

where the subscript x may be v, c, or a, and jS indicates the type of time integration, as discussed above,

is modeled using equations (25a,b,c). Note that = = 0 and Ps = Pv = P© at t = 0.

Some material is converted to gases during pyrolysis; the rate as which gases are created is given by

equation (25e). The gases are lost from the cell, and they carry away all the energy they contain, i.e., there is an enthalpy loss which does not affect the temperature of the remaining mass. Moreover, since

the gas diffusion equations have not been explicitly included, any possible loss of heat from the escaping

hot gases to the cooler solid in other regions of the substrate is ignored. Thus, the rate of heat gain in

the cell due to pyrolysis is (almost exactly) given by

S^ = LV{R,H,^R^^H^^-R^^H^^-i>^cT) Watts (64)

The Rj are given in Section II. B. 3, and dpj/dt is given by equation (25e).

7. Time Integration

A choice must be made for the time integration method. Following the discussion of Belytschko (1983,

pp. 55, 419, 445), the advantages of explicit time integration are:

(1) fewer calculations per time step;

(2) simple algorithm logic and structure are simple, implying that it is good for testing new

ideas;

(3) ease of handling complex nonlinearities;

(4) requirement of little core storage compared to implicit methods using direct elimination

procedures; and

(5) high reliability in terms of accuracy and completing the computation.

The only notable disadvantage is that explicit time integration is only conditionally stable so that a very

large number of time steps may be required.

With regard to accuracy, since implicit methods are unconditionally stable, they can easily be used with

too large a time step, leading to significant time integration errors. The stability requirements for explicit

time integration force the time step to be so small that the time integration error is almost always smaller

than the spatial discretization error. Of course, it is also possible to use a spatial discretization that is much too large.

The addition of pyrolysis to the model required significantly smaller cells in the region of interest, and

therefore required significantly longer execution time. Hence a better method than explicit time

integration was required. The following method attributed to Saul'yev as described by Larkin (1964) and

Clausing (1969) was adopted:

24
Again consider the one-dimensional presentation of Figure 1. Assume that the calculation of cell temperatures is proceeding in the positive x direction. Then at cell i, Tj.j is a known quantity and

can be used in computing T; n+j. Equation (51) becomes

At/, , = C, ,T. ,-C, T, = LUq, ,+q, +S. .) (65a)

During the next time step, calculate cell temperatures in the negative x direction. In that case equation

(51) becomes

At/, , = C, ,T ,-C, T, = Lt(g, +Q, , + .) (65b)

Both equations (65a) and (65b) are unconditionally stable because of the inclusion of the implicit terms; operating together, the truncation errors during successive time steps tend to cancel leading to an O(At^)

algorithm. This method is directly expanded into 2 or 3 dimensions by adding the j and k position

indices and the qj., qj^, q^,.,, and q^.+ heat gains. The key factor is that it is only necessary to solve

implicitly for one cell temperature at a time. There are no time-consuming simultaneous equations to be

solved. Tests indicate that this method is very accurate except for the possibility of some small

oscillations as with the Crank-Nicholson method.

A question remains on when to evaluate Sj, Cj, and /Cj (which is implicit in the q's) in equation (65a).

Numerical errors are minimized by evaluating Sj at time step n-f-V^. These terms are all functions of

temperature in the SUBSTRAT program. Furthermore, they are such complicated functions that the

equations cannot be solved directly. The following equation must be solved implicitly for Tj n+ji

* ^ij.k-.n*l ^ij.k*.n'^ ^i.J.k.n*p )

where

^U.k.M ~ ^iJ,k,m^^i,J,k^i.j.k.mTi,i,k,m

with m = n and n-l- 1. The total solid density in the cell is the sum of the virgin, char, and ash densities

[equation (25d)] and is a function of time due to the pyrolysis reactions:

The cell pyrolysis for the time step is evaluated at T = (T^ -t- T„+i) / 2, and component densities also

at jS = Vi. Equation (66) was rewritten in a form allowing numerical solution by a standard secant method (Conte, 1972, pp. 299-306). This completes the outline of Larkin's method. The net result of

using Larkin's method is that the time step is no longer limited by grid size, and pyrolysis is modeled

with a reasonable execution time.

25
D. EXPERIMENTS

1. Experimental Arrangement

The next step is to compare the results predicted by SUBSTRAT with experimental results, in which an

electrical heater was substituted for a cigarette. The schematic of an experiment designed to measure time

to ignition is shown in Figure 4. This consists of the heater element from an automobile cigarette lighter,

fitted with a concentric jacket to permit an air purge. The purpose of the air purge is to keep evolved

products from the sample from being ignited to flaming, rather than to merely smolder. The heating

element is raised to varying temperatures, all of them high enough for the element to be glowing (500 - 900 °C). The purging jet comes through the jacket, past the face of the element, and then out normal

to it, toward the sample. It picks up a good deal of heat as it travels through the device. The resulting

flux distribution beneath the heating element is shown in Figure 5. As we would expect, it is axially

symmetric; it is well fitted by a Gaussian profile. Note that the higher the temperature the disk is heated

to, the higher is the peak flux.

The measured flux to the gauge is not the same as what the substrate sees, however, because the latter heats up, whereas the flux gauge does not; call the former <f>^ and the latter, 0^. It follows from equation

(3) that

= <l>.<^*4>. = 4>^*M7;-r,) (68)

where Tg is the temperature of the hot purging gas, T^, that of the (cold) gauge, and <^s » <t>^J^. For the case where the substrate is being heated, however, the convective contribution goes down:

= 4>ra^*'»[7;-7;(0] (69)

where Ts(t) is the temperature of the substrate, which increases continuously. Thus the total flux decreases monotonically. Although we do not know the gas temperature Tg, we do not need it; for, we

can combine equations (68) and (69) to obtain

4),(0 = i>,-h[T^it)-T^] (70)

The heat transfer coefficient h was discussed in Section II. B. It is found for this experimental

configuration in Appendix D. Analysis of the experimental results, given in Appendix D, yields a number of values important for this experiment, including h and the disk temperatures; the results are shown in Table D-1; h is found to be a weak function of 0exf

When a mock-up consisting of flexible PU foam covered by #12 cotton duck was placed in position under

the heat source, smoldering ignition occurred after a certain amount of time; the ignition delay depends on the intensity of the flux, as indicated qualitatively by equation (19), Section II. B. 2. The experimentally obtained ignition delays are plotted in Figure 6 as a function of the peak heat flux.

26
Heat Source

Gas Purge

Figure 4. Schematic of heat source for ignition tests

27
FRONT/BACK HEAT FLUX SCAN

15 MM LEFT/RIGHT; 5.4 MM VERTICAL

1.0 ' ' i > 1 10.0 10.5 11.0 11.5 12.0

FRONT/BACK DISTANCE (CM)

Figure 5. Flux profile from heat source, as measured by a total heat flux gauge

28
The ignition delay times as a function of the external flux are given in Table 2.

Table 2. Substrate Ignition Delay Times

</>,,,(kW/m2) 18 25 34.5 44

ti«(s) 472 70 37 22

Since <f>J^g^ is constant for each run, then according to equation (19), a plot of (^g^t versus t-t'^ should

yield a straight line whose intercept is the critical flux, (f}^^^^. When this is done, however, me result is not a straight line. More generally, then, we assume that tfie relationship has the form

(71)

Only the correct choice for <|)^^^^ will yield a straight line in a logarithmic plot of vs tjg. When that

is obtained, the slope of the line is -p; it was found that 0j.rit — ^^-^ kW/m^ and that p = 1.087, with

A varying from 757 to 816; the average value is 798. If we take p = 1, A increases slowly with 0^^^^.

The fact that p « 1, rather than 1/2, indicates that the substrate behaves like a thermally thin material;

that, in turn, suggests diat it is principally the fabric which is involved in the heating and ignition. If we

assume that the principal cooling mechanism is radiation, then <|)^^^^ = 16.9 kW/m^ implies that if the

radiative absorption coefficient of the fabric is a = 0.9, then = 759 K « 485 °C ~ a remarkably

high value. Even with a = 1.0, Tjg = 739 K = 466 °C, considerably greater than the 400 °C measured

independently. This shows that there is substantial convective cooling during the heating/ignition process.

2. Material Data

In order to see whether the program can calculate tjg correctly, it is necessary to calculate the surface

temperature under the heater. In order to do that, however, it is necessary to have correct input data. Among these data are the thermophysical data, Kpc.

Fabric . First, consider the fabric. The fabric that was used in the experiment was #12 cotton duck. Material data for cotton, especially as a function of temperature, are surprisingly difficult to find, even though it is a common and long-used material.

In Appendix E, we find that reasonable values for the thermophysical constants for this particular cotton

fabric (#12 cotton duck) are

p = 620 kgjm^

and

29
(S) 3IAII1

Figure 6. Time required to ignite the fabric, for different (initial) heat flux exposures. Crosses

correspond to the discussion in the text. Filled circles correspond to a different set of

experiments.

30
K(T) = 0.28505 K,(r) + 0.84554 K^„(r) W/m-K (72)

where

K,(T) = x^(TJ (T/TJ = 1.242x10-^7 W/m-K (73)

(with T in Kelvins) and k^^^ =1.7 k^^^, with

1 +

The constants are: a = 2.6464x10"^, b = 245.4, and c = 12. Then k(25 °C) = 0.1435 W/m-K and c(T) = it(T)/pa « 7819/c(T); thus, c « 1122 J/kg K at 300 K.

Char. The cotton decomposes and pyrolyses to form char. According to Parker (1985, 1988), the

specific heat of char is just about that of carbon:

Cc(T) = 1.43 + 3.55xl0-^T - 7.32x10'*/t2 J/g-K, with T in K (75)

The thermal dilfusivity of wood char is approximately constant: « 2.1x10"^ m^/s (76)

The density of the char depends on whether the fibers contract while pyrolyzing, or not. Finally, the

thermal conductivity of char is found from

ic,(T) = aePcCcCT). (77)

Reaction Kinetic Parameters . When the analysis of a pyrolyzing material is carried out, it is done for a thin layer, in order to ensure uniformity of temperature and of oxygen concentration throughout the

sample. It is then simplest to find the reaction rate in the form

expi-EJRT) min"^ (78)

for the degradation reaction, rather than as in equation (20); here W^ = W(i(t) is the (instantaneous)

weight of the remaining (virgin) material = weight of sample minus weight of char and ash, and W^ is the original weight. Note that is commonly given in reactions/minute. The actual reaction rate is given by

The relationship between (20) and (78) is simple: since = WjA^j and = WqA^q, then as long as VoA^d = const, WjAVq = Pd/po, and, as in equations (20),

p;exp(-...)

where

31
(in the appropriate units). For pyrolysis in pure nitrogen, Ohlemiller (1991) found: = 7.49x10^'* min-^ n = 0.60, and = 182.4 kJ/mol

Instead of these parameters, however, we shall use the results of Kashiwagi and Nambu (1992) for a

cellulosic paper pyrolyzed in air; that will maintain internal consistency in the reaction set. The global

kinetic constants given by Kashiwagi and Nambu for the thermal degradation reaction of their cellulosic

paper are

Ed = 220 kJ/mole

Ad = 1.2x10'^ min-*

n, = 1.8

Hd = - 570 J/g

The parameters for the other two reactions are given in their Tables 1 and 2, reproduced below. For the

oxidative pyrolysis reaction,

\ 0/

'/CP exp

KT

(79)

where X^. is the volume fraction of oxygen. The kinetic constants are:

E = 160 kJ/mol = 1.5x10'^ min-i

Hop = 0.5

nfop= 1-3

and: H„p = 5.7 kJ/g

For the char oxidation reaction.

with the kinetic constants

C C OX exp

[ m

(80)

E^ = 160 kJ/mol

A. = 3.4x10" min-'

Hco = 0.78

He = 1.0

and = 25 kJ/g

The fraction of ash which remains is 9% by mass. Equations (79) and (80) transform to the standard

form in the same way that (78) does. Thus

Kp = ^op^o '^"Pa

32
Heats of Combustion. Brandrup and Immergut (1989) give the heat of combustion from several

measurements; as expected, the results vary:

Hc(cotton) « 18850 J/g

H<,(fabricl) « 15450 J/g

Hc(fabric2) « 16700 ± 250 J/g (fabric weight was 180 g/m^)

Note: as material pyrolyzes, we should take c(T) and k(T) for the combination of materials; for the sake

of simplicity, however, we assume that the virgin fuel, the pyrolyzed material, the char, and even the

ash, all have the same specific heat. Then the appropriate density to use in equation (12), for example,

is the density of the solid, p^, which is the sum of the bulk densities of all the components of the solid

(see equation (25d)). As pyrolysis and combustion proceed, decreases monotonically, assuming no

fabric shrinkage.

Foam. For the padding component of the substrate, we have PU foam, for which we have the following

data (Ohlemiller, 1991):

p = 0.032 g/cm^ = 32 kg/m (81a)

Cp = 1.46 J/g K (81b)

and /e(T) = 0.03613 + 2.003x10-'* T, (81c)

with T in "C.

There is one other datum for PU foam: <Kpc> , the mean value of the thermal inertia, was measured

in the LIFT apparatus (Quintiere et al., 1983, Quintiere and Harkleroad, 1985, Quintiere, 1988). However, the foam melts before it ignites, so that p (and, consequently, k) go up substantially; moreover, the

measurement entails all the other model approximations. Hence the measurement is not directly useful.

Foam Kinetic Parameters . The first pyrolysis reaction, which results in the collapse of the foam

structure, yields 0.7 g of liquid and 0.3 g of vapor for every gram of foam which pyrolyzes. The

reaction rate, RRj a (mass gasified/min)/(mass of foam), is given by

RR^ = 0.30 B

m ( exp (82)

where m = 1.5 B = 3.36x10^'', and

Eb = 187 kJ/mol.

Also, Po = starting density,

p = current density, and

Pr = density of residue « 0.70po.

These densities are those of the solid phase {p^ of order Ig/cm^), rather than that of the foam (p^ « 0.03

g/cm^). The total heats of pyrolysis of PU foams range between 400 and 800 J/g; for this foam, it is closer to the upper value. Assuming it varies linearly with the mass loss, we can estimate 30% of 800,

or

33
Hp ~ 240 J/g

Finally, for the diffusion coefficient for O2 in the fabric, see equation (28).

E. RESULTS AND DISCUSSION

Rather than give the final results immediately, it was decided to describe the evolution of our thinking,

including some of the false steps we took, because a good deal can be learned that way, and it could be

instructive for anyone else who might work on this problem in the future. With that in mind, we first consider the results of calculations made assuming an inert substrate.

1. Inert Substrate

First, the computer program was used to calculate the peak surface temperature under the heater as a function of time for several peak fluxes. The results are shown in Figure 7,

The ignition temperature of cotton was measured to be 390-400 °C (Ohlemiller, 1991), According to Figure 7, if Tjg « 400 °C, then ignition is attained for the 44 kW/m^ exposure at t « 18 s, and for the

34.5 kW/m^ exposure at t « 50 s. These values compare very favorably with experimental values of

22 and 37 s, respectively.

On the other hand, ignition is not attained at all for the 25 and 18 kW/m^ cases, while the calculation

predicts the substrate will reach the ignition temperature, albeit at long exposure times. A conceivable reason might be that the ignition temperature given above was overestimated. In order to intersect the

<^ = 18 kW/m^ curve at t « 472 s, the ignition temperature would have to have been about 290 °C; it

is very unlikely that such a large error in the measurement of T^^ would have been made. Even if it were, that would yield ignition times of 7, 12, and 25 s, respectively, much shorter than the measured

times of 22, 37, and 70 s.

[Three possible reasons why the calculated ignition times might be so short are:

(a) endothermic pyrolysis taking place, which slows the temperature rise;

(b) Significant radiative heat transfer within the material, so that the flux is absorbed in depth, rather than mainly at the surface; and/or

(c) Tjg = 290 °C and the estimates of /cpc(T) were in error, the actual value being larger.

Now, if we consider endothermic pyrolysis, then we must also consider exothermic pyrolysis. Second,

the eflFect of radiative heat transfer will be small when T is near T^, and small in comparison with the

efl'ects of pyrolysis when T> >Ta. Finally, consider item (c): If Kpc were larger, then the rate of rise

of Tg would be smaller, as can readily be seen from equations (17)-(19). The final temperature, however,

is independent of Kpc. Hence not only would Tj^ have to be much smaller, but the estimates of Kpc

would have to have been 9 times too small. This is all possible, but extremely unlikely.]

34
2. Results with Pyrolysis

On the other hand, we know that pyrolysis must take place. Figure 8 shows the effect of including

pyrolysis for the Q = 25 kW/m^ case. Curve a indicates what the temperature of the fabric surface

would be if the fabric (and foam) were inert. Curve b is the result of "turning on" the (endothermic)

degradation reaction. Note that although the rate of growth of temperature is slowed down quite

perceptibly, as expected, the temperature the surface reaches asymptotically is exactly the same (assuming

the same surface absorptivity/emissivity as before).

However, the exothermicity of the oxidative pyrolysis and of char oxidation are far greater than the

endothermicity of thermal decomposition, and should overwhelm it. We cannot, therefore, consider

thermal decomposition alone: if we consider pyrolysis, we must include all the major reactions. Another

way of seeing this is to note that in order to have ignition, we must have a "thermal runaway," where

exothermic pyrolysis heats up the material faster than heat diffusion and surface losses can carry the heat

away.

Curve c indicates what takes place when the oxidative pyrolysis reaction is included, as well: it is exothermic, and adds a great deal more energy than the degradation reaction removes. The result is, logically, that the temperature rises rapidly, almost like a thermal run-away. However, the temperature

reaches a peak, than declines. The reason is simple: the fuel is rapidly exhausted. Once that happens,

the heat source is reduced to the original external flux. This is clearly depicted in Figure 9, which shows

the fuel density as a function of time. The ordinate is temperature, in °C, and density in kg/m^. The

curves marked Tq and pq correspond to the cell with the highest temperature. Those marked T| and pj

correspond to the (laterally) adjacent cells, and those marked T2 and P2 to the next ring of cells. These

calculations were carried out assuming no char oxidation. We see that as the reaction accelerates, T^

"runs away" and the density plummets from its virgin value to that of char. Most significantly, the peak

temperatures develop just after the density falls.

Finally, curved of Figure 8 shows the result of adding the char oxidation reaction as well: that additional

heat source takes the pyrolysis "over the top." The temperature continues to run away. (We have

arbitrarily cut off the calculation at 600 °C, here.) The question then arose: why was there an oscillation

in Tg? The qualitative explanation is that if a cell is too large, then the surface/volume ratio is small, and

the heat cannot diffuse away rapidly enough. This is confirmed by Figure 10. With a cell size taken to be a bit less than half as large, the amplitude of the oscillation declined considerably, and with the cell size halved again, the oscillations have almost disappeared.

Presumably, as we continue to decrease the cell size, the numerical errors should become progressively

smaller, until a further decrease in cell size would make no further difference in the results. However,

it was found that the results were apparently not converging with decreasing Ax, and that the sensitivity

varied with Q. Investigation of the reason for this great sensitivity revealed that it lay in our use of a

variable grid size: whereas for grids of constant spacing, the numerical approximations are correct to second order in Ax, that accuracy drops to somewhere between first and second order. Indeed, if the

numerical error is proportional to (Ax)", then n is given by a complex formula (Torrance, 1985, pp. 43

and 51). The "flavor" of that expression is given by n ~ (2+a)/(l + a), where a » dAx/dx, TTiat is, the accuracy depends inversely on the rate of growth of the grid size. An explicit calculation is given

in Appendix A; see Table A-1 there.

35
Figure 7. Peak surface temperatures of substrate as a function of time, for the four exposures

36
700

100-

t(s)

Figure 8. Peak surface temperatures of substrate exposed to Q = 25 kW/m^, as a function of t, for

a. Material assumed to be inert d. Char oxidation also included, as well

b. Thermal degradation only e. Only one reaction: oxidative pyrolysis

c. Oxidative pyrolysis as well

37
800

"^0 20 4-0 60 SO 100 120

t(s)

Figure 9. Temperature and density of central cell (subscript 0), adjacent cells (subscript 1), and cells in next ring around center (subscript 2), as functions of time

38
Figure 10. Temperature of central surface cell (i.e., peak temperature) for three different grid sizes: 1.25,

0.50, and 0.25 mm cubes. Char oxidation was purposely left out

39
soo

700

600-

500-

^0 10 20 30 40 50 60 70 SO 90 1 0C

t(s)

Figure 11. Peak temperature for the 25 and 34.5 kW/m- cases, assuming (a) no pyroiysis, and (b) all three pyrolytic reactions

40
By reducing the rate of increase of grid size, therefore, the accuracy was increased, and the results made

to converge better. Moreover, it was decided to switch from an explicit solution method to Larkin's

semi-implicit method (see Section 11. C). After those changes, we arrived at the results shown. The

results of four calculations are shown in Figure 11. These were made

with an assumed heat transfer coefficient h = 20, and an assumed oxygen concentration of 20%. For

initial fluxes with peak values of 25 and 34.5 kW/m^, calculations are first made assuming that the

substrate is inert, and then pyrolysis is "turned on." We see that for the 34 kW case, the asymptotic

temperature lies at about 430 °C, and that a thermal runaway begins at about 300 °C, at t » 20, and is completed at t « 25 s; that would dien be the ignition time. It is difficult to state precisely what "the

ignition temperature" is, in this case (but see below).

For the 25 kW case, the asymptotic temperature lies at about 360 °C. The first perceptible deviation of

the curve for the reactive case from that for the inert material occurs at T « 280 °C, and is clearly

established at 3(X) °C. The thermal runaway takes place at t = 48 s. These times may be compared widi

the experimentally-established ignition times, which were about 22 and 70 s, respectively. Thus we have

achieved semi-quantitative agreement. In Figure 12, on the otiier hand, we see that a thermal runaway,

and therefore ignition apparently did not take place for the 18 kW case. Here, the asymptotic

temperature is about 290 °C; this is apparently not high enough to produce a runaway.

Wherein lies the difficulty? Most likely, one or more of the input values is incorrect. Altering jcpc

would principally change the time scale. Figure 13 shows the effect of changing the assumed oxygen

concentration: the highest curve reproduces the upper curve in Figure 12. The next two curves show

what happens when <Y> is assumed to be 0.15 and 0.11, respectively. Finally, the bottom curve

(again) corresponds to the inert case. Thus, it is not [O2] having been chosen too low that prevents

ignition, and one or more of the kinetic parameters is probably in error. In Figure 14, curves a and b

again reproduce Figure 12 (on a different scale). For curve c, the char oxidation rate was doubled. It

is apparent that the curves overlap completely. A preliminary conclusion inferred from this was that the

observed result was due to all the char that is produced already being oxidized. Observation of the

kinetic constants for oxidative pyrolysis and char oxidation, however, makes it clear that the latter is three

orders of magnitude slower than the former. Therefore merely doubling the char-oxidation rate will only

perturb the energy output slightly ~ so slightly that it will not even show up in the figure.

For curve d, the oxidative pyrolysis rate was doubled (the pre-exponential factor A was doubled),

doubling the char-production rate tfirough this branch; this indeed produced a thermal runaway. Curves

e and f are the results of increasing A by 20% and 10%, respectively. Thus, a quite modest increase in A predicts ignition, although at 230 s, rather that the measured 472 s. Such an increment is not only well

within experimental error, but is entirely plausible, when the likely differences between the cellulosic

paper and a cotton fabric (with different impurities) are considered. On the other hand, it was assumed

that y = 0.20. For y = 0.15, A must be increased by 30% in order to get ignition (the resulting T(t) curve is very similar to the "best" one). Although this is greater than the 10% increase found above, it

is still entirely plausible.

We have so far considered the sensitivity of the results to the oxygen concentration, die thermophysical

constants of the fabric, and its kinetic parameters. Surprisingly, there are two other significant

parameters: first, if we use h = 22 for the heat transfer coefficient in equation (69) (as suggested by the

results of Appendix D, shown in Table D-1), rather than the assumed h = 20, the "asymptotic"

temperature (that at t = 500 s) reaches only 272 °C, rather than 292 °C, and there is not the faintest possibility of achieving ignition, unless is substantially smaller than 160 kJ/g. One solution is to assume that h = 20 is the correct value to use, since the theoretical calculations in Appendix D could

41
350

50-

"^0 50 100 150 200 250 300 350 400 450 5O0

t(s)

Figure 12. Peak temperature for the 18 kW/m- case, with the same assumptions as in Fig. 11

42
43
800

Figure 14. Peak temperature for the 18 kW/m- case, with various assumptions for the pyrolysis:

Curve a, no pyrolysis; curve b, "standard" pyrolysis; curve c, double the char oxidation rate; Curve d, double the oxidative pyrolysis rate; Curve e, 1.2 times the oxidative pyrolysis rate; Curve f, 1 . 1 times the oxidative pyrolysis rate

44
easily be off by 10% or more. Another resolution is possible, too: the second parameter which is important in this threshold region ("threshold," because 18 kW/m^ is close to the critical flux, 16.9

kW/m^) is the thermal conductivity of the foam padding. In all the runs made above, it was assumed that

K = 0.056 J/m K and c = 1.9 J/g K, for the foam at T = 20 °C. If it were assumed that k = 0.096,

instead, then the foam would act as a more efficient heat sink, and the surface temperature could be

expected to drop; indeed, a calculation showed that the peak surface temperature at t = 500 s fell to 272 °C for the inert fabric. In fact, nowever, the value 0.056 for the thermal conductivity was for a foam of density 48 kg/m^! Transforming that for a 32 kg/m^ foam, according to equation (62), yields

0.036, almost exactly what is given in the present test (see equation (81c)). This reduces the heat sink,

and must yield an increased asymptotic surface temperature. Using the foam parameters given by

equations (81), and h = 20, the t = 500 temperature indeed rises, from 292 °C to 300 °C; with h =

22, it is 294 °C, and we only need to increase A by 10% to get runaway.

The "best" set of parameters for the fabric and foam, then, is given in Section II.D.2. For the heat

transfer coefficient, use the values in Table D-1. With that set, we obtain the curves shown in Figure

15. The corresponding calculated ignition times are given in Table 3. The calculated values are all about

half the measured values, and the trend is well-predicted.

Note that the polyurethane foam begins to melt and recede from the fabric when its temperature reaches

about 300 °C, thereby decreasing the heat-sink effect of the padding, and accelerating the heating of the

fabric. This effect has not been included in the model.

Table 3. Comparison of Substrate Ignition Delay Times

Ignition delay, tj^ (s) Q (kW/m^) Calculated Measured

18 209 472

25 37 70

34 19 37

44 13 22

It has been suggested that we might avoid the necessity of explicitly including the pyrolysis reactions by

choosing some effective ignition temperature. This is not, in fact, feasible: if we take the measured

ignition times and mark them on the four curves in Figure 7, corresponding to the inert assumption, we

find that they intersect these curves at widely varying temperatures; see Table 4. It is apparent that this "ignition temperature" is a strong function of the external flux.

45
Figure 15. Peak temperature as a function of time, for all four cases, using the best set of input data

46
Table 4. Surface Temperatures Which Would be Attained by the Substrate

at the Measured Ignition Times If the Substrate Were Inert

4>,^ (kW/m^ ) 18 25 35

1 1 44

291 342 380 400

Ti,(K) 564 615 653 673

On the other hand, if we define "the ignition temperature" as the point on the Tg(t) curve where the

temperature is rising at some rapid rate, e.g., 100 °C/s, then from Figures 11 and 13 we see that that gradient is attained at the approximate temperatures shown here:

Table 5. Surface Temperatures at a 100°C/s Rate of Rise

(Data From Figures 11 and 13)

0,^, (kW/m2) 18 25 35

Ti« (°C) 390 390 380

Ti.(K) 663 663 653

These calculated temperatures are close to each other, and close to the 390-400 °C which has been

measured.

F. SUMMARY FOR SECTION II We have created several computer programs, of which the central one is SUBSTRAT2. This program

calculates the temperature history at every point in a substrate (the upholstered furniture) which is subjected to a strongly localized heating flux on its top surface (assumed to be horizontal). The (solid)

substrate consists of two layers, the top one being a fabric and the lower a foam pad; there may be a thin

intervening air gap. The substrate is taken to be a rectangular parallelepiped, and it is broken up into

several thousand cells. The ambient oxygen level can be set at whatever value one wishes.

The solid is subjected to a nonuniform and time-varying heating flux at its top surface, and (simultaneously) experiences convective and radiative heat losses. Moreover, the solid can undergo pyrolitic reactions;

we consider three, here: one endothermic step (thermal degradation to char), one oxidative pyrolysis to

char, and oxidation of the char (to ash). The way T(r,t) is found is by solving the PDE which describes

the diffusion of heat in a solid, equation (12), numerically. The equation set is very stiff, because of the

highly nonlinear form of the (Arrhenius) reactions. We have therefore used a semi-implicit method to solve the equation set.

If and when the temperature rate of rise at some location suddenly "accelerates" to a value high enough

that the surface glows (that is, T > 500 °C or so), we can say that smoldering ignition has occurred.

The program will not tell us whether flaming ignition takes place. It also does not treat the case where

the flux is applied in a crevice, such as is formed between the seat cushion and the seat back. The

47
program does not take into account the effect(s) of the foam possibly melting and receding away from

the fabric. Finally, it also does not take oxygen diffusion within the cushion explicitly into account; hence in certain threshold situations, where a small change in oxygen concentration determines whether

ignition does or does not take place, the results are ambiguous. What is and is not in the program is listed in Table 1, section II.B.

It is important to note that it is often difficult to obtain the needed kinetic and/or thermophysical

parameters for the material; or, when available, to know how accurately they are known. Therefore this caveat must also be made: even if the program were perfect, its results are only as good as the input

parameters which are supplied. On the other hand, it should accurately reproduce (or predict) trends.

There is a user-friendly front end for the input, described in Section II. C. The program runs well on a 386-Ievel computer with a math co-processor, and faster on a 486-chip computer.

An experiment was carried out to ignite the substrate. In trying to reproduce those experimental results,

it was found that the calculated results are sensitive to the input values chosen, especially the kinetic

parameters. It was found that the preexponential factor found by Kashiwagi and Nambu (1992) for the

oxidative pyrolysis reaction in a cellulosic paper had to be increased by 10% for the cotton duck fabric

in order to get ignition for the lowest of a set of heating fluxes to which the substrate was exposed. The

result was semi-quantitative agreement with the observed ignition times.

48
ra. CIGARET, A MODEL OF A QUIETLY SMOLDERING, ISOLATED

CIGARETTE

A. GENERAL DESCRIPTION

A cigarette consists of small strands of cured tobacco leaf (varying in length from about 0.5 to 12 mm,

with the average about 3-4 mm), held in place by a paper wrapping. The shape is usually cylindrical;

the radius of the cross-section is R. The circumference, 27rR, generally lies between 17 and 25 mm.

There is usually a filter at one end.

A schematic illustration of a smoldering cigarette in its quiescent phase is given in Figure 16. The section

marked C in this figure is char, most of which is oxidizing at a rate sufficient to make it glow. The peak

temperature in this region is some 800 to 850 °C; that is, about 1 100 ± 25 K. The glowing coal is shaped approximately like a thick cone. The cone length varies; it is usually on the order of 1 cm.

The heat created by the char oxidation diffuses, convects via hot gases, and radiates away in all directions. Part of this heat is carried back towards the virgin tobacco, the zone marked VT in the figure.

The tobacco just behind the burning coal then dehydrates and decomposes; this is the region marked P

in the figure (P for "pyrolyzing"). This is the region from which a visible plume of smoke rises. At the

front of the cigarette is the residual ash, marked A in the figure. EA stands for "evanescent ash"; that

is the ash which generally falls off" a cigarette, or is flicked off by the smoker. Since we will eventually assume that the cigarette rests on a substrate, the part EA might not, in fact, disappear.

We will see in Section III.C that a good deal of detail is lost when the model is assumed to be onedimensional. Therefore the present model is two-dimensional; i.e., cylindrical symmetry is assumed.

This is not a bad assumption for the cigarette freely smoldering in air, but is not so good when the

cigarette lies on a substrate. We nevertheless assume that whatever asymmetries are produced by that interaction can be simulated by spreading the efl'ects of the substrate out evenly over the circumference.

B. DYNAMICS OF A SMOLDERING CIGARETTE

A good deal of efl'ort, both experimental and theoretical, has gone into understanding the (isolated)

smoldering cigarette. In this Section, a brief, qualitative discussion of the dynamical processes which

occur will be presented and the salient results from earlier experiments noted. It is not intended to be

an exhaustive review. For a more detailed description of cigarettes, see Gann et al. (1988).

1. Pyrolysis Rates

There have been many careful studies of the pyrolysis of tobacco; we will only consider some global

results here. A detailed description of cigarettes can be found, for example, in Section 2 of Gann et al. From Table 2-3 there, it is found that for those experimental cigarettes, the average mass, exclusive of

the filter, is 870 ± 470 mg. The variation for commercial cigarettes is lower. The mean mass loss rate

is 60 ± 20 mg/min, corresponding, on average, to some 14 min of smolder time. This agrees with the

common observation that an unattended cigarette will smolder for 10 to 20 min.

49
Figure 16. Schematic of a smoldering cigarette. The zone marked VT

contains the virgin tobacco. P is the pyrolyzing region; the

section marked C is char, A is the residual ash, and EA stands

for "evanescent ash," explained in the text.

8mm

DISTANCE FROM LINE OF PAPER BURN (mm)

Figure 17. Isotherms in the gaseous part of the cigarette during quiet, steady

burning. The shaded region indicates peak reaction rate; i.e.,

the glowing coal.

50
Sandusky (1976) shows that if we start with one gram of tobacco, then after pyrolysis about 400 mg of

charred residue remains. This char includes 150 mg of inert material, which remains as ash upon

combustion. That is, about 70% (600 mg) of the combustible portion (850 mg) of the tobacco evaporates

and/or pyrolyses away, leaving about 250 mg of material which can burn completely. To avoid

confusion, we will refer to that portion as "char," and the char plus ash (total — 400 mg) as "charred

tobacco residue," or CTR, hereafter. We will not go into any more detail here regarding pyrolysis rates.

We now try to determine a sensible approximation to the stoichiometry of the reactions. This will be

useful when the production of CO, H2O, and CO2 are modeled. The contents of the tobacco are quite

complex, as shown in Sandusky (1976) and elsewhere. The reactions are complex as well. A simplified

summary of the reactions is given at the bottom of page V in Sandusky. TTiis is here simplified still further, by assuming that the dehydration and pyrolysis processes leading to char have exothermic as well

as endothermic terms which cancel exactly.

When the char bums, the principal products are CO and CO2, leaving about 150/400 = 37.5% ash, by

weight. This residue is somewhat higher than the amount left by the tobacco analyzed by Baker (1975).

Sandusky's data is used here, complemented where necessary by Baker's data. Char is assumed to consist of a mixture of carbon and of one other compound which will yield C/H/0 ratios comparable to those given by Baker, and for which we know the heat of combustion, H^.. A model molecule is mucic

acid, CgHioOg. This has the molar heat of combustion = 2021 kJ/g mol (Weast, 1976). Its molecular weight is 210; thus h^. = 9620 J/g. The atomic percentages of char (CTR minus ash) are given by

Baker as:

The "model" char is then

char = xC + yC6Hio08,

where x and y are such that the above ratios are satisfied. From the carbon and hydrogen fractions, x « 12 and y « 0.5 are obtained; nitrogen is ignored. Upon burning the char,

12c + 0.5 CgHioOg + n02 -» (15 - f) CO2 + f CO + 2.5 H2O

In order to determine f, one more relationship is needed. This is obtained by examining the energy

produced. We assume complete combustion of the mucic acid,

0.5 CgHioOg + 4.5 O2 3CO2 + 2.5H2O,

and incomplete combustion of the carbon:

12c + (n - 4.5)02 -* (12 - f)C02 + f CO

The heat of combustion of C (to CO2) is 393.5 kJ/g mole (32.77 kJ/g), while the heat of combustion of

CO is 283.0 kJ/g mol (23.57 kJ/g), so that the heat of partial combustion (C to CO) is the difference,

110.5 kJ/g mol. With these assumptions, the entire char produces

2. Stoichiometry

0 = 25.33%

C = 69.40%

H = 2.18%

N = 3.09%

51
G = 12(393.5) - 283.0 f + 0.5(2021) J

The molecular weight of this model char is 249. Thus the heat produced in combustion is

h^. = G/249 J/g

The mean heat of combustion of the tobacco, as measured by Sandusky, is

<h^> = 17360 ± 210 J/g.

Therefore, f = 5.0 ± 0.18, and n = 11.75.

Thus the CO/CO2 ratio is f/(15-f) = 0.5. This is typical for smoldering, though very large in

comparison with the ratio in flaming combustion.

Note that even though there is only 1 mole of CgHigOg for every 24 moles of C, its heat of formation

contributes 23% to the overall heat of combustion, a quite significant fraction of the energy balance. The

overall model reaction, finally, is

12c + 0.5 CgHioOg + 11.75 O2 - IOCO2 + 5C0 + 2.5 HjO (83)

3. Regression (Burning) Rate

This study only considers the quiescent phase of smoldering, i.e. , between puff's, since that is the situation

when the cigarette lies on the substrate. The paper in commercial cigarettes is chemically treated so that

it decomposes and burns at the same velocity as the tobacco during these quiescent periods; this occurs

in a region on the order of 1 mm in width. This location is called the "paper burn line," and will serve

as the origin of coordinates (x = 0) in our further work. It is assumed here that the reactions in the

paper do not influence the tobacco reactions, although we realize that there may well be both thermal and

oxidative interaction.

Cigarettes vary within a factor of about two in the speed of propagation of the smolder wave. The

average cigarette tobacco colunm is about 7.3 cm long, and it takes 15 ± 5 minutes to be consumed,

without drawing or puffing on it. The burning rates have been measured to lie between 40 and 85

mg/min, consistent with the value given in Section III.B.l. The regression velocity is thus

= L/t = 5.4 ± 1.8 mm/min.

The variations depend, at least, on the cigarette radius, the packing density, whether the tobacco leaf has

been expanded or not, the moisture content of the tobacco, and the kind of paper, its permeability, and

how it has been chemically treated. Note that the moisture content does no! affect its peak temperatures

(M. Samfield, 1986). Rather, it affects the burning velocity, because it takes time and energy to evaporate the water.

52
4. Reaction Rates

As will be seen in Section in.D.2, the reaction rates in combustion are well represented by Arrhenius

expressions, which can give enormous rates at high temperatures. What limits the reaction rates is primarily the availability of oxygen; that, in turn, depends upon the rate at which the oxygen can be

transported into the combustion (Ohlemiller, 1985). Therefore, the rate at which oxygen is allowed to diffuse into the cigarette is the principal determiner of the burning rate, unless the coal is so cooled

locally as to cause limitations from the kinetics. The diffusion of oxygen will be discussed further in Section in.B.9.

The mean reaction rate for the cigarette is

where 0.85m is the combustible fraction of the tobacco mass, V is the volume of the tobacco column,

and T the total duration of the smolder. We use a nominal mass of m = 1 g. The volume is V = ttR^L,

where the average L = 7.3 cm and R « 0.4 cm; finally, t » 840 s. Then <RR> « 0.28 mg/cm^s.

If the reaction rate is described by an Arrhenius expression, then it must be highly peaked in the hightemperature region. Indeed, the simple approximation of an infinitely high reaction rate occurring over

a surface which is a paraboloid (or conoid) of revolution might be expected to have some validity. This

approach was used by Gugan (1966). Since the actual reaction rate is finite, this "surface" is in fact a

thin region . This shell-like region is shown in Figure 17 (taken from Baker (1975)) by the shaded

paraboloid. The highest temperatures occur at die centroid of this shell. Thus, suppose an observer were

positioned at x = 0, i.e., at the paper bum line. As time progresses and the smolder wave moves to the

left, the peak reactions occur in a contracting ring, starting at the periphery and contracting to a point,

leaving ash on the outside of the ring.

We have seen that the mean mass loss rate is about 60 mg/min for a 870 mg cigarette; hence something

like 70 mg/min for a 1 g cigarette. Of this, about 70% is evaporated/pyrolyzed with low energy production or loss, 30% is lost by char oxidation, a highly exothermic reaction. That is, about 20 mg/min are being oxidized this way. With R^o ~ 17 kJ/g, we end up with an energy production rate of about 340

J/min, or about 5.7 W. Note that with the variations among cigarettes described above, this value could

easily be 50% greater or smaller:

<RR> = 0.85m/VT,

5. Heat Production Rate

<Q> « 5.7 ± 2.8 W

6. Energy Balance

This heat input is balanced by losses; without going into detail, we can estimate:

a. Enthalpy loss via outgassing: 2.0 ± 0.4 W

b. Convective losses from the surface: 1.7 ± 0.4 W

c. Radiative losses from the surface: Qr « 2.1 ± 0.4 W

53
As we see, although the energy losses balance the energy production, the uncertainties in these estimates are very large indeed.

d. Conductive losses: These apply when the cigarette is in contact with a substrate, and are discussed in Section IV. For the "free" cigarette (that is, the cigarette quietly smoldering in wind-free

air, distant from any substrate or wall), there are no conductive losses.

7. Distributions Within the Cigarette

Temperature. A representative temperature distribution in a smoldering cigarette is shown by the

isotherms in Figure 17. Figure 21, also taken from Baker (1975), shows the measured volume percentage

of oxygen. Note that this value is essentially zero in the region of maximum temperature, so that the

reaction (oxidation) rate there is low, in spite of the high temperature.

In order to be able to calculate the heat flux delivered to the substrate, the distribution of surface

temperature must be known. In particular, most important is the peak surface temperature, Tp.

Unfortunately, it is not an easy quantity to measure. Indeed, it is not even a well-defined quantity, since

(as is seen from Figures 17 and 19) the gas and the solid temperatures are not quite the same.

Moreover, the result depends on the measurement technique. Baker's measurements yielded peak surface

temperatures of about 550 °C. Egerton et al. (1963) found 616 °C for Tp. Lendvay and Laszlo (1974),

using an IR technique, found that Tp = 600 °C. We will take this to be the mean value. Since peak

surface temperatures may well vary by 50 °C among cigarettes, the peak surface temperature is 600 ±

50 °C.

Oxygen . The oxygen concentrations in the freely smoldering {i.e., away from the substrate, in quiescent phase) cigarette are shown in Figure 18. It can be seen that the smoldering is indeed

cylindrically symmetric. In fact, in our model it will maintain that synmietry even when interacting with

the substrate. Consistent with an oxygen-diffusion-controlled process, the oxygen concentration in the center drops to very low values. The peak reaction rates generally occur in the region where the mole

fraction of oxygen, x(02), is less than one percent. One can write down some plausible analytic expressions for the radial distribution; but the actual distribution yields some surprises, as can be seen from

Figure 18. For example, although one would a priori expect the concentration to increase monotonically

from the center towards the boundary, examination of the actual dependence along a slice at z = -1-6 nmi

shows unexpected maxima and minima.

Gas Velocity . Some simple estimates show that the radial outflow velocity of hot gases from the

cigarette can be as high as 7 cm/s. Thus inward oxygen diffiision is a counterflow problem. Indeed, this suggests that the principal inflows would be through the ash, away from the paper-burn line, and then

axially towards the reaction zone. Note that the buoyancy of the hot gases will in fact affect the gas flow.

54
Surface Heat Flux . Consider Figure 19, giving the measured temperature of the cigarette surface,

parallel to its axis. The net convective flux from this surface to the ambient is

4>^(x,t) = hmx,t) - TJ (84)

Similarly, the net radiative flux is

4>,(x,0 = ecO[T\x,t) - r/] (85)

where is the emissivity of the cigarette surface, a is the Stefan-Boltzmann constant, and is the

ambient temperature. These fluxes, when integrated over the cigarette surface, must yield the energy loss

rates Qr given in Section III.B.6.

The heat flux from the escaping hot gases, i.e. , the enthalpy loss, can be calculated once we have Up the

radial convective velocity, as a function of x. Upon integration over the surface, it then yields the total enthalpy loss, Qe8. Pressure

Upon drawing on a cigarette, a person will generate a (negative) pressure corresponding to several

hundred Pa (several centimeters of water, where one atmosphere corresponds to 101 kPa and 10.4 meters

of water.) In the quiet state, the pressure diff'erence is positive, but a minute fraction of this value. This

pressure difference can be estimated:

The flow through a packed bed follows Darcy's law,

m" = fLp (86)

But

(87)

where Pg is the gas density at the surface and u^ is the radial velocity. Pg is given by the ideal gas law.

Hence,

Pg = pJ,/Tg (88)

U Q T

An = -Ll^ (89)

With u^^ « 0.07 m/s, p^ = 1177 kg/m^ and Tg « 900 K « 3T^, we have:

. 7x1.2x10-2 o , Ap « Pascals

3/

Sandusky (1976) estimates that for peripheral flow, f « 10"^ s/cm. Hence:

Ap « 30 Pascals

One atmosphere is 101 kPa, so this corresponds to 3 x 10"^ atm, or about 3 mm of water equivalent.

55
Oxygen

+ 12

Figure 18. Oxygen concentration in the freely smoldering cigarette during

quiet, steady burning.

DISTANCE FROM LINE OF PAPER BURN (mm)

Figure 19. Isotherms in the solid part of the cigarette during quiet, steady

burning.

56
9. Gas Transport

There is a boundary layer surrounding the cigarette, within which the oxygen concentration drops as the

surface is approached. Fresh air is supplied by the ambient via convection and diffusion, as indicated

by the arrows in Figure 20. That the oxygen must transport inward is clear, since the volume of air needed to consume the cigarette is about 1000 times the volume of the solid. Thus, 99.9% of the oxygen

needed for combustion of the cigarette must come from the ambient. Some part may be drawn in from

(through) the substrate, when the cigarette is in contact with it, but that is estimated to be a small fraction

of what is drawn in from the atmosphere Ohlemiller et al. (1993, Appendix C), and it has not been

modeled here. If there were no substrate, the inward diffusion would be approximately symmetric, the

only deviation from symmetry being caused by the upward gas flow due to buoyancy.

Although the cigarette paper is permeable, it largely inhibits the diflFusion of oxygen from the air into the

tobacco, so long as it is intact. Once it has burned, its resistance to oxygen diffusion is much diminished.

Since no measurements of this residual resistance have been published, we assume that it is zero. Therefore, most of the oxygen which is needed to sustain combustion diffuses into the cigarette in front

of the paper bum line. [Not all of the needed oxygen is drawn in that way, however. If the paper is made non-porous, the cigarette will self-extinguish.]

Similarly, the considerable gradient in CO2 which is established by combustion means that there is simultaneous outward diffusion of CO2. Stoichiometrically, each molecule of O2 reacts with a single

carbon atom. This is how the cigarette loses most of its remaining mass upon combustion (recall that 60% of its mass, at a given point, already disappeared during dehydration and pyrolysis). Some of this mass is convected out, and some of it diffuses out, as just described. It is not difficult to show that the

convection rate is one to two orders of magnitude greater than the diffusion rate. Note, to begin with,

that before there is any counterflow diffusion of oxygen, about 60% of the tobacco column mass

evaporates/pyrolyses away. This process occurs just adjacent to the combustion region. Thus at least 2/3 of the total outflow is due to convection.

10. Diffusion Coefficients

The diflFusion that we are concerned with is that of O2 in N2. We will call the diffusivity of O2 in N2,

Dq. For porous media, the effective diflFusivity depends on the porosity, or the void space. It is shown

in Szekely et al. (1976) that

^^ffl^c = 0.677$^-^* (90)

for f» < 0.7. Although the total void fraction is <i> = 0.85 (Muramatsu, 1981), we use this relationship

nevertheless. We thus obtain the effective diffusion coefficient

Deff « 0.56 Do = 0.112 cm^/s

at ambient temperature. This is precisely the value used by Muramatsu (1981). The diffusion coefficient

is assumed to be the same for all the gases.

57
Substrate

Figure 20. Schematic of how air enters the smoldering cigarette

when it is resting on a surface.

0

\

^

\

\

\

\

\

Instant of ^

Contact ^ <

1 V,

A

B

Figure 2 1 . Schematic of the energy losses of a smoldering cigarette as a

function of time, when it is dropped onto a surface. Curve A

corresponds to a dense, inert substrate. Cuve B corresponds

to a substrate which undergoes exothermic reactions.

58
From equation (16.3-1) of Bird et al. (1960), we find that for O2 in N2, the temperature dependence is (ji \1.823 —-— cmVs (91)

293.16

j

We shall use the same coefficient, but the theoretical temperature dependence

= Do(T/273)^^^ (92)

rather than the empirical relation above.

11. Conduction Losses

These apply when the cigarette is in contact with a substrate, the situation of concern in actual fire initiation. According to the estimates made in Gann et al. (1988), the initial power loss to the substrate

is about 2.1 ± 0.35 W, going down to about one-third of this in the steady state.

C. PREVIOUS MODELING EFFORTS

A number of attempts have been made to model the isolated smoldering cigarette or similar object. All

of these attempts make some simplifying assumptions in order to make the problem tractable.

An early, related model is due to Moussa et al. (1977). They experimented with, and then modelled,

the smoldering and extinction of a cellulose cylinder without any paper wrapping. They assumed the

smoldering to be steady-state, and they treated the problem as one-dimensional {j.e., no radial gradients).

They found, experimentally,

• the smoldering velocity to be closely related to the maximum temperature in the cylinder;

• reasonable agreement between theory and experiment for the extinguishment limit;

• the rate-limiting step in the combustion to be diffusion of oxygen to the char, which is in good accord with experiment.

They also calculated values for v, the propagation velocity of the smolder wave. However, that calculation depends on an uncertain parameter, and the overall accuracy of the model is questionable.

A much more realistic and detailed model of a cigarette was produced by Sandusky (1976) and by

Summerfield et al. (1978). It is also a one-dimensional model, but it considers the steady-draw case, rather than the free-bum condition. The model is heterogeneous; that is, it explicitly takes into account

the fact that the tobacco comes in shreds. It considers a two-step process: pyrolysis to a char, followed

by the oxidation of the char. It ignores water evaporation. The char-oxidation reaction is assumed to be a linear function of the oxygen concentration. The combustion model is time-dependent rather than

steady-state, and considers the cover paper indirectly, via a varying surface permeability to oxygen. It uses a sophisticated treatment of heat transfer inside the cigarette, including heat transfer by radiation as well as by solid phase conduction.

59
The model consists of ten simultaneous, coupled partial differential equations (PDEs), and the starting

condition assumes the presence of fixed amounts of ash, char, and tobacco. The calculation predicts the

burning velocity fairly well, as a function of draw rate. The dependence on oxygen mole fraction in the

atmosphere is less well predicted, although this is not important for the current purpose. The calculated

gas temperature profile is fairly good: the peak is about right, but the width of the distribution is too narrow. Since it is a steady-draw model, diffusion and natural convection within the cigarette are considered as being of relatively minor importance.

The most elaborate model is that due to Muramatsu et al. (1979, 1981); it was developed in two stages. In the 1979 reference, they developed a one-dimensional, steady-state model for the pyrolysis of the

cigarette. They focused on the evaporation-pyrolysis zone in a naturally smoldering cigarette. The model

considers:

• pyrolysis of tobacco obeying Arrhenius kinetics,

• evaporation of water from tobacco following a mass-transfer and rate-determined

process,

• weight loss of tobacco due to pyrolysis and evaporation,

• internal heat transfer characterized by an eff'ective thermal conductivity which includes

approximate radiation heat transfer,

• heat loss attributable to free convection and radiation from the outer surface of the

cigarette and to endothermicity of the evaporation process, and

• smoldering speed.

These processes are expressed by a set of simultaneous ordinary differential equations that are solved

numerically by the Runge-Kutta-Gill method. The equations are ODE's rather than PDE's, since they

are independent of t and depend only on x, the position along the axis. The propagation velocity, v, is imposed . They do not include the convection or diffusion of gases other than water vapor. They take

the existence of the paper wrapping into account only through its eff'ect on the loss of water vapor.

These approximations work well. The model yields good agreement between theory and experiment for

the temperature and density along the axis in the pyrolysis-evaporation region. Thus, the agreement for

T(x) and p(x) is good for x < 0, i.e., before the char-oxidation region, which is not considered in this part of their model. [As in Section IILB, we use x = 0 as the boundary between the pyrolysis region

and the char-oxidation region, at the surface.] For x>0, the calculated temperature profiles deviate

substantially from measured ones, as might be expected. The dependence of the profiles on the imposed

velocity shows only semi-quantitative agreement.

In the 1981 Muramatsu work, a char-oxidation model was added to describe the processes occurring in the region x > 0. The model is quite detailed; it takes two char oxidation reactions into account and is two-dimensional (cylindrically symmetric). Unlike the Sandusky model, it is a homogeneous model, i.e. ,

it does not take directly into account the fact that there are solid particles and a gaseous medium. Energy

loss is through radiation and convection at the outside surface of the cigarette. Again unlike the Sandusky

model, it assumes that there is no temperature diff'erence between the solid and gaseous phases. This is not a bad approximation for natural smolder, and simplifies the problem considerably. Heat transport

by thermal radiation inside the cigarette is taken into account in a somewhat diff'erent way than is done

60
in Summerfield et al. The thermal conductivity at any point is assumed to be isotropic. Similarly, the

temperature-dependence of the gaseous diffusivity is taken into account explicitly. Finally, the

pyrolysis/evaporation model and the char oxidation model are tied together through an energy-flow

matching condition at the pyrolysis/char-oxidation boundary to obtain the appropriate smolder velocity.

The results of calculations made by Muramatsu for six representative cigarettes agree well with

experimental data:

• The peak temperatures, when expressed in °C . are only 2.7% to 5.7% higher than the

experimental data. [Since the temperature calculations are made with T^ « 20 °C as the

datum, it is appropriate to use degrees Celsius for comparisons.]

• The smolder velocities are 14% high, on the average, varying between 4% low and 26%

high.

• The calculated variations of smolder rate (V) and peak temperature (T^) with R (the cigarette radius), die packing density, Pp, and the moisture content in the tobacco shreds are

indistinguishable from the experimentally observed variations. The dependence of V and

T^J^ on ambient oxygen partial pressure is not well-predicted; this is not an important

consideration in our study.

• The calculated distributions of temperature and oxygen concentration in the char

oxidation region are in agreement with measurement, except for a scale factor: the

predicted distribution is narrower than observed. Because of this last point, Muramatsu

et al. 's model cannot be used directly to obtain the flux emitted to the substrate: the too- narrow temperature distribution would substantially underpredict the energy output of the

cigarette to the substrate. However, die model is excellent for some purposes; in particular, to estimate cigarette smolder velocities.

Although the model is very good, it may be appropriate to list some of its limitations:

• Prior to decomposition, the paper wrapping is assumed to be impervious to oxygen.

According to experiment, this should result in the cigarette going out.

• Perhaps in order to be consistent with the above assumption, radial convection of gases

is not included.

• Their calculation uses an iterative procedure to converge to a solution. If the calculation

has not converged after 1000 iterations, the unconverged result is nevertheless accepted

as correct.

Mitler (in Gann et al., 1988) and Mitler and Davis (1987) developed a detailed, homogeneous model of

a ft-eely smoldering cigarette, called CIG25. Unlike Muramatsu's models, it is time-dependent rather than

steady-state and two- rather than one-dimensional. Unlike Sandusky's model, difl\ision and natural

convection within die cigarette are included, since there is no steady-draw convection to overwhelm these

eff'ects. We may make some further comparisons; CIG25 diff'ers from the model introduced by Sandusky

et al. in a number of ways:

61
• Sandusky's model is heterogeneous, with tobacco shreds embedded in a gas. Therefore

the local "ambient" for oxygen varies, and the rate at which oxygen reaches the shred

surface is given by

^n(yo-ys) (93)

where y^ is the local ambient mass fraction of oxygen and y^ is the value at the shred

surface.

CIG25 is homogeneous, with the cigarette core being a mixture of "solid" and "gas."

The oxygen available to react with the solid is the local value, obtained by solving the

diffusion equation for y, the oxygen mass fraction.

• The Sandusky et al. model is a steady-draw model with substantial axial convection and no natural smolder; CIG25 considers natural smolder only, and neglects axial convection

in comparison to radial convection.

• The Sandusky model is a one-dimensional model which averages over the cross-section.

CIG25 is two-dimensional, having radial gradients.

• CIG25 takes only one reaction into account, effectively considering a char cigarette,

rather than a real tobacco cigarette.

Of course the models have a number of similarities; for example, they both neglect the production,

transport, and effect(s) of water.

Mitler (1988) also described a "semi-empirical" model. This simple version consists of a number of

correlations, partly based on results of making parametric runs with Muramatsu's model and corrected

by experimental data where possible.

D. MODELING THE CIGARETTE

1. Assumptions

Based on the successes of the prior models, this Section lists the assumptions and equations valuable to a good, yet tractable model of a cigarette on a substrate. CIGARET uses a subset of them and is described in Section III.D.3. The following are facets to be included and their physical forms:

1. The cigarette model is two-dimensional, with axial and radial coordinates. It is also timedependent, yielding V directly, whether or not it is constant. If a steady state were

assumed, then the equations would simplify; however, it would then be necessary to choose the smolder velocity correctly [i.e., as an eigenvalue), as Muramatsu did. Moreover, if the equations are taken to be time-dependent, the heat and mass transfer

equations are parabolic partial differential equations (PDE's) which are easier to solve,

in some ways. Finally, the convergence for the steady state equations is very slow.

Indeed, Muramatsu found that one thousand iterations would not suffice, at times.

2. The cigarette is modeled with one pyrolysis reaction and one char-oxidation reaction.

62
3. The (internal) gas flow velocities depend on the pressure gradients; since the axial gradients, of order Ap/L, are generally much smaller than the radial gradients, of order Ap/R,

axial convection is neglected, as being much smaller than radial convection.

4. The water (pre-existing or produced during combustion) in the tobacco column is also

ignored.

5. The tobacco column is treated as a homogeneous, uniform mixture of gas and solid; there are no tobacco shreds.

6. The gas and solid phases are at the same temperature, locally.

7. Species or temperature gradients widiin the tobacco shreds will be neglected, consistent

with assumption #5.

8. Radiation transfer within the cigarette is treated as an efl'ective thermal conductivity.

9. The thermal conductivity is the same function of temperature throughout the cigarette,

i.e., whether in ash, char or tobacco.

10. The paper behavior appears only in the boundary conditions.

11. Consistent with assumption #10, any (axial or radial) gradients in the paper are ignored; an eff'ective mass transfer coefficient is used to model diffusion through the paper,

12. The gases generated or heated by combustion move radially outward to the side boundary

(aside from diff\ision). Thus there is a radial flow calculated strictly by mass conservation. The gas pressure within the cigarette is assumed to be only negligibly different

from atmospheric, and therefore no momentum equations are written.

13. The gas phase is quasi-steady; that is, dp^/dt = 0,

14. No particulate aerosols are produced.

15. The cigarette combustion zone retains its cylindrical symmetry when lying on a substrate. This is observed to be approximately correct for some cigarettes on some fabric/foam

substrates, especially after the cigarette "recovers" from the transient cooling eff'ects of

the substrate.

16. Prior to its ignition, the substrate's presence has two eff'ects on the cigarette, both of

which can be assumed to apply symmetrically (on the average) to the entire cigarette

periphery:

(a) The oxygen supply to the cigarette is reduced by some factor. The degree of

reduction will depend on the permeability of the substrate.

(b) The thermal eff'ects {e.g. , as a heat sink) can be calculated in a decoupled manner

(this assumption is weak, and may eventually have to be dropped).

63
It should here be made explicit that this model is for an isolated (that is, not in contact with any surface),

quietly smoldering cigarette. Such a model will predict, among other things,

• the external heat flux from it, and the extent of the heating zone,

• the velocity of smolder propagation, and

• how these depend on: the radius of the cigarette (R), the tightness of packing (via the

void fraction <j>), the type of tobacco {via its thermophysical and kinetic parameters: heat

of gasification Hy, heat of combustion H^,, activation energies and pre-exponential factors,

density p, thermal conductivity k, specific heat Cp, etc), and the wrapping paper, via its permeability and its ignition, or decomposition, temperature. A better model would

explicitly include paper thickness, chemical composition, porosity, and kinetic parameters, as well.

2. Governing Equations

The equations which describe the mass and energy transport in a smoldering cigarette will now be

presented. Any simplifying assumptions beyond the 16 above will be indicated as each equation is discussed.

It is necessary to clarify some of the terms appearing in the equations. In the decomposition of the

tobacco, one gram of tobacco, upon pyrolysis, produces n,. grams of char. The other 1-n,. grams are gaseous products. Each gram of char reacts with no2 grams of oxygen to produce some heat, n^ grams

of ash, and 1 - n^^ + nQ2 grams of gaseous products of combustion. Since the combustion process mostly

proceeds at very low oxygen concentrations, it can be expected to be incomplete, and nQ2 is smaller than

the stoichiometric value. In the following equations, x is the axial coordinate.

Mass Conservation.

dp, (94) = -Kl -nc)R^ ^ (1 - n^)Rc^] .

where is the mass density of the solid, Rp is the pyrolysis rate (in gm/cm^s) and R^^ the char oxidation

rate. Axial convection is dropped, according to assumption #2, and the equation of continuity in cylindrical coordinates becomes

|-[p,(l - *) + P,*] + -^(<l>rp,«,) = 0 (95)

dt * r dr '

where r is the radial coordinate, Pg is the mass density of gases, <f> is the void fraction in the cigarette

(i.e., the volume fraction of gas, rather than of tobacco shreds), and u^ is the radial velocity of the gas. Since the shreds are assumed not to shrink, <i> remains constant. Assumption #1 1 implies

and using equation (95), the gas equation results:

64
Momentum. The procedure used in CIGARET is to assume no pressure difference so that the ideal

gas law is satisfied. Mass conservation then gives the velocities directly, via the gas equation, Equation

(97). The reason is that any attempt to calculate the convective outflows and velocities from a calculated

pressure difference would require that we know the coefficient in Darcy's law (Section III.B.9) very well,

and that we accurately calculate very small differences of large numbers. That would be extremely

demanding of the computation and possibly prohibitive.

Species . V is the total volume of the cigarette; the total volume occupied by the solids in the

cigarette is

7, = (1-4>)K

We define the mean densities

p. =m./K, i=A,C,T (98)

where A, C, and T stand for ash, char, and tobacco, respectively. In this equation, m^ = total mass

of ash, etc. Then the cumulative density of the solid parts of the cigarette is

Ps = Pa + Pc + Pj (99)

Each of these densities varies with time.

The equations for the changes in tobacco and char densities are

= -R (100)

dt

and

dpc

C p CO

The ash is inert and simply accumulates, so that an equation for is not required.

Oxygen. We must know the oxygen concerntration everywhere in order to calculate the charoxidation rate. The other gases are inert or nearly so. Moreover, the gases and solids are assumed to have the same temperature, so there is no energy exchange, and we do not need to consider the CO2

diffusing out. Therefore, of the gaseous species equations, it is only necessary to include that for oxygen.

Rather than dealing with the oxygen density Po2» equations are best written in terms of y, the oxygen

mass fraction:

y - Pojpg (102)

65
The equation is

^*ar ^' 'dr dx[^' 'dxj rdr[ ^' ' dr)

(103)

where is the oxygen diffusion coefficient, and the axial convection term has been dropped.

Energy . Assuming that the gas is quasi-stationary, as in assumption #12 above, and that the

specific heats of all gases are the same and independent of temperature, then the energy conservation

equation can be written in terms of the temperature. It is:

(1 -<^)p.C, ^jws s — + <|)p„ M,C„ — = —\k— + \rk — ^yg r g dx\ dx) r dr\ Br

(104)

where = specific heat of the solid

Cg = specific heat of the gases

k = thermal conductivity of the cigarette (see assumption #8)

Q^o = energy released from char oxidation (lower heat of combustion)

Qp = energy absorbed in (endothermic) pyrolysis

The internal heat transfer has a radiative and a conductive component, as described earlier. The

expression used here is the same as that used in Muramatsu (1981) and is due to Kunii (1961); for porous

materials, ,j

k(T) = (1 - *2/3)jt^ ^ $1/3^

g 3

» p

.XV

^ (105)

where h, is a heat transfer coefficient for radiation:

h^ = €^4oT' (106)

Dp is the mean pore diameter, $ is the total void fraction (including the void space in the shreds, and is therefore larger than <^), Cj is the emissivity of the shreds, kg is the thermal conductivity of the gas, and

kg that of the solid shred (and depends on the shred's mass density).

Since the gas pressure in the quiescent cigarette is very nearly the ambient air pressure, and since there

is not a great difference between the molecular weight of the product gases and air, the ideal gas law

permits one to write the gas density in the form

f>g = PgoLIT (107)

where T is the absolute temperature.

66
Reactions . Finally, expressions for the tobacco and char reaction rates are needed. It is assumed that each is given by an Arrhenius relation:

= p;Z^exp(-£;^//?D (108)

where the exponent m is to be determined experimentally, as is the "frequency factor" (or "preexponential factor") Zp. Ep is the activation energy for the pyrolytic reaction and R is the universal gas

constant. The tobacco density pj may be expressed as

Pt = Ps y? (109)

where y^ is the tobacco mass fraction.

Similarly, the char-oxidation reaction rate is taken to be

Rco = Pc p4 ^co ^^V(-EcoIRT) (110)

where (again) the exponents n and p are to be determined experimentally, and the densities are written

in terms of the respective mass fractions:

Pc = P8yc and Po2 = Pgy (111)

Initial Conditions . Before ignition of the cigarette, the initial conditions are

Ps(x,r,0) = Pso = Pto (112a)

Pc(x,r,0) = PA(x,r,0) = 0 (112b)

T(x,r,0) = To (112c)

y(x,r,0) =ya(-yambient = 0.232) (112d)

For the calculations, however, it was assumed that a match had been applied to the x = 0 end of the

cigarette, producing some reactions. These used up only a small fraction of the tobacco, but most of the

oxygen that had been in situ. The initial temperature distribution was assumed to be very high (1000 K)

in the first millimeter of the tip, then to decrease linearly to ambient temperature in the next millimeter.

The oxygen mass fraction behaves in complementary fashion, as follows. For the following calculation,

we may ignore the stoichiometry developed earlier and simply assume that the combustion of char

proceeds according to

C + [(1 *a)l2\0^ aCO^ + (1 -a)CO

Then 16(1 + a)/ 12 grams of oxygen are needed to bum one gram of carbon (char). As combustion lowers

the local mass fraction of oxygen from y^ to y, it will lower the relative density of char (and hence

tobacco) from 1 to x, where

X = 1 - [0.75/(1 +fl)n,](p,/Pr)(y« - y)

This permits us to enter the initial density of tobacco so that it is consistent with the assumed distribution

of y(02) y. Po is the actual mean tobacco density in the virgin cigarette.

67
Boundary Conditions . On the axis,

dy]

dr

=0 foraUxyt (113)

/r=0

At the lit end of the cigarette, the temperature boundary condition is:

k = €^a(r^ - T^) + /r,(r - (114)

for all r and t. T = T(0,r,t) and = TambientHere, too, h^ = convective heat transfer coefficient (at the end), a is the Stefan-Boltzmann constant, and = emissivity of the cigarette at the x = 0 end. We are assuming a grey body. The value to be used

here should be that of the tobacco ash- It is assumed that the ash never falls off the cigarette, not realistic for active smoking, but appropriate for the case when the cigarette rests on a substrate. Moreover, if it were noj made, the geometry would be continually changing, and the boundary conditions would become

exceedingly complicated.

There is usually a filter at the other (x = L) end. It has been observed that the presence of a filter increases ignition propensity (Gann et al. 1988). This is probably due its limiting axial flow of air through the tobacco colunm. This model presumes no such flow. Thus, for the temperature boundary

condition at the other end.

(115)

for all r and t (where now T = T(L,r,t)).

Next, the cigarette's side surface must be considered. It is covered by paper. This paper wrapper in principle should also be included with its own set of equations, which include its reaction kinetics.

However, the amount of heat released when it burns is negligible in comparison to that released by the

tobacco. If a unique paper decomposition or ignition temperature can also be specified, then the paper's

kinetic equations can be replaced by two appropriate boundary conditions, which simplifies the problem

considerably.

If there were no paper, then the temperature boundary condition at the sides would be

- kix,R) ^-^j

= €,(x)a(r* - r/) + h{T - (116)

The emissivity of the cylinder surface is difl'erent for ash, char, or virgin material (except at the tip, the

relevant material is paper, rather than tobacco, as is pointed out below); thus the emissivity is a function

of X. Also, the thermal conductivity at the surface will in general depend on whether it is ash, char, or

virgin material. The assumption which is made here is that k is the same for all three (assumption #9).

Similarly, although the heat transfer coefficient at the sides, h, will be difl'erent from what is at the ends,

we simplify by assuming that

he = h

68
Moreover, since the ends generally emit much less heat than do the sides, it is not so important that we

use the exactly correct values there. The emissivity has one value for the virgin region, another for

the char region, and a third for the ash region. Note that these refer to paper char and ash, and may

occur at different temperatures, and therefore different locations, than for the just-underlying tobacco char

and ash. Because of these differing emissivities, c^. will vary along the surface. Although ejix) may be

a continuous function, the simplest approximation that can be made is that the paper pyrolyses, ignites,

and disappears at some paper ignition temperature Tjp, so that there are just two values for e^. - that for

the virgin paper and that of ash. These meet at the paper burn "line." Experiments (Sandusky, 1976)

have shown that Tjp « 450 ± 100 °C.

The peak surface temperature occurs towards the "tip" end of the cigarette (i.e., at x < 0), and therefore

it is most important that be accurately chosen.

Finally, the oxygen boundary conditions must be considered. It is assumed, for the sake of simplicity,

that the filter prevents any oxygen diffusion at the cold end (x = L):

D >1 =0 (117)

where D is the diffusion coefficient for oxygen within the cigarette. At the other end, the conditions must

be

(118)

where y^ = y(0,r,t) is the oxygen mass fraction at the "hot" end of the cigarette and y^ is the ambient

fraction, defined in equation (112d). D is a function of temperature; for the sake of simplicity, it is assumed that it is the same function for ash as it is for virgin tobacco and for char. 7 is the mass transfer

coefficient, sometimes referred to as kg or k^; see, for example, equation (93). Muramatsu (1981) gives

it as

= 6.38 X 10

-3 T^.75(^T- T)(T+ 123.6)

1/4 cni/& (119)

for the effect of the boundary layer, where T is the mean value between T^, (the local cigarette surface

temperature) and T^ (the ambient temperature), and all temperatures are in Kelvins. This holds for the

free cigarette, where there is no blockage by a substrate

For the side surface, the counterflow convection must be taken into account:

dr

y(ya - y,)

x<x.

(120)

where Xp = Xp(t) is the position of the paper burn line at time t, u^ is the (outward) radial velocity at the

surface, and y^ is the oxygen mass fraction at the surface.

The paper resistance to oxygen diffusion can be expressed as

69
7p = Dp/6 (121)

where Dp is the diffusion coefficient in the paper and 6 is the mean paper thickness. The total resistance

of both paper and boundary layer is then

7 = (7b"' + 7p-')-' (122)

According to equation (120), the remains of the paper wrapper, where the paper has burned, present no

barrier at all, and the radial oxygen diffusion at the boundary depends only on the properties of the

boundary layer.

It is instructive to make estimates for the magnitudes of these terms. From equation (106), we see that

at the peak surface temperature, 600 °C, Dg « 0.856 cm^/s. For the gradient, we see that x(0^ « 0

to about halfway out along the radius, so that we can take

^ > Zil^ (123)

dr RJ2

For 7, we use y « 0.8 cm/s (from Table 5-B-l, Gann et ai, (1988)). Then equation (120) relates u, and y^. From Figure 18 (Figure 5-3b in the same reference), [OJs = 8 ± 2% in the most active region.

Therefore, y^ « 0.09 ± 0.02. Hence,

2D - 1 0.856

0.2

- 0.8

( 0.232 _ A

\ .09 ± .02 J

- 1 = 2.9 ± 0.5 cm/s

If the gradient is indeed given by equation (123), then in order to be able to satisfy equation (120), we

must have

Ur < 2D/R. (124)

That is, too large a radial convective flow will prevent any oxygen from being able to enter the cigarette

through the sides. This quantitatively confirms our physical intuition.

3. Choices for CIG25

The model CIG25 used a subset of the above equations which was believed to capture the most important

processes. The key simplifications were:

• Since the pyrolysis is only weakly endothermic, it was dropped ~ that is, it was assumed

that Rp = 0 in equations (97), (101), (103), and (104).

• Probably more limiting, it was assumed for simplicity in some parts of the program that

Pj. and Pg are constant and uniform. In other parts of the program, the ideal gas law was

properly used; hence, the relationships were inconsistent.

• There was no boundary condition for convective flow. This was not a serious issue,

although the intact paper forms a barrier against convection. However, since there is no

pyrolysis in this model, the only place where gas flow can originate is where high

70
reaction rates occur, producing high temperatures and expanding gases. All this takes

place largely in front of the paper burn line, and therefore one would not expect signifi- cant flow behind that line. Therefore there was no need for such a boundary condition.

• The radiative losses from the surface were linearized (see equation (131)).

Because this was only a "partial" model, we would expect the results obtained with the correct values of

input parameters to be somewhat unrealistic. The input parameters therefore have to be modified, in

order to compensate for this.

E. NUMERICS

1. Discretization of the Equations

For the numerical calculations, the cigarette is divided into 10 or 20 cylindrical shells, i.e., Ar = 0.4 or

0.2 mm, and into slices of the same thickness. Ax = Ar. This last equality makes the expressions a little simpler and more compact.

Originally, the equations were discretized using central diff'erences for the spatial derivatives. This gives

simple expressions which are accurate to second order in Ax. However, it is easy to show (e.g., Peyret

and Taylor (1983), Chapter 2) that this approximation introduces an artificial (i.e., numerical) diffusion

term with a negative sign. Hence, if the actual diff\ision term is not large enough, numerical oscillations

will begin, and instability ensue.

The stability criteria for the homogeneous convection-diflFusion equation

^ ^ - DV'f=0 (125)

dt dx dy

are

and

Lt <. — (126)

\A\' IBP

Af^Al!, (127)

42)

assuming that Ay = Ax. Equation (126) is the Courant criterion.

Care was taken to satisfy the constraints on At and on the magnitude of the mass diffiision coefficient. Nevertheless, severe oscillations arose within a few hundred time steps. This was traced to the

inhomogeneous source terms. Thus, when the convective-diffusive equation we are solving has a source

term, that invalidates the criteria for At found above. The source terms are of Arrhenius type, and

produce considerable stiffiiess in the equations. It would have been necessary to constrain At to a microsecond or so in order to avoid these oscillations. These problems could have been overcome by

using implicit solution schemes and/or operator-splitting methods (Wichman, 1991). Another approach

which is often used is to quasi-linearize the equations. This, however, would have reduced the accuracy

71
of the solutions. Thus, it was decided to use the sHme numericEl technique as is used in CIG25. It is described in Section III.E.2.

If one takes Ax = 0.2 mm, then with the present DIMENSION statements, that only allows for cigarettes

of 30 mm length (ordinary cigarettes average about 73 mm). However, it is likely that a quasi-steady

state or ignition would develop long before 30 mm of cigarette was consumed.

2. Method of Solution

Ames (1969) shows that an equation of the form

* f(xj,u)^ + g(x,t,u) = p(x,t,u)^ (128)

can be put into the Crank-Nicolson form as follows:

-^6l[ui*' + u/] + -j-f[ih,(j^^l2)k,(u{*' + «/)/2]6,(ur* + «/) +

g[ih,(J*m)k,(ui*' + tt/)/2] = p[i/i.(/+i/2)*.(«r' + «/)/2] ' * '

(129)

where 6^ is the central-difference operator, h = Ax, and k s At. The equations above are precisely of

this form. The method used to converge to a solution at each time step, is to iterate according to an

under-relaxed Gauss-Seidel scheme. See Mitler (1988) for more details.

F. Improvements in CIGARET Over CIG25

There are five categories of improvements in CIGARET:

• The physics has been improved.

• The input is much more user-ftiendly. It also accepts input from SUBSTRAT.

• The output is in a form which can be used directly by SUBSTRAT.

• The program has been redesigned to be quasi-interactive with SUBSTRAT and can

calculate the effects a substrate has on the smoldering.

• The documentation is improved.

These are described in the next three subsections. There is one more improvement which does not fit into any of the above five categories: The power of personal computers has increased substantially in the

past six years, so that it is now feasible to run the program on a 486-level computer. Therefore the

72
program has been modified to run on a PC, rather than on a CYBER, as was the case for the antecedent

program, CIG25.

1. Physics

The improvements in physics are four:

1. The gas density is now given by the ideal gas law, rather than taking it to be constant in

the oxygen and the gas velocity equations, as was done earlier. See equations (97) and (103). In CIG25 the approximation was made that the product PgD^ could be taken from inside the

derivatives in the first two terms on the right-hand side of equation (103). This eliminated the

effects of gradients in the equations, and also allowed Pg to be factored out from all terms in that equation but the last. That is, in the identity —V-(p^DVy) = DV^y + rVyvf^^^] (130)

the second term on the right-hand side (where the ideal gas law was used to eliminate Pg), was

dropped. This term has now been inserted into the program. Note that here we have taken into

account not only the gradients in the gas density, but those of the diffusion coefficient as well.

2. Likewise, the char density now behaves correctly in the velocity subroutine, rather than

assuming also to be constant. In equation (97), the reaction rate R^q depends on the char

density (which falls with time). In CIG25, the density was taken to be constant, in this equation.

Similarly, so was the gas density, Pg. Both of those densities now vary with time and position,

properly.

3. The view factors for radiation exchanges have been calculated (see Section IV) and

incorporated into CIGARET.

4. The radiation losses from the surface of the cigarette are now correctly calculated. That

is, they are as given by equations (114), (116), and (148), rather than using the linearized form

Kss ' h'(T^ - T^) (131)

where h' is an effective heat transfer coefficient, which only approximates the effect of radiation

loss.

2. Input and Output

The input section has been redesigned for much greater ease of comprehension and use. For example,

the program now requests dimensional physical quantities, rather than the ratios used in the program and

called for in the input section of CIG25. The program has also been redesigned so that it can indirectly

interact with SUBSTRAT. See the USERS' GUIDE (Section II.H), for further discussion of these points.

The output files have been changed to accomplish the following:

73
• There is a complete "dump" at the end of each output time (file CIG0UT2), the previous

dump being eliminated. This permits one to resume a calculation which has been

interrupted for whatever reason. It also keeps the file from being too large.

• A small subset of CIG0UT2, the axial and the surface temperature distributions, is saved

at each time step, and added to the other output file, CIGOUTl. SUBSTRAT uses only

the surface distribution part.

It is appropriate at this point to indicate the modifications which have been made to TMPSUB2 to produce

SUBSTRAT:

• SUBSTRAT uses the cigarette surface temperature distribution computed by CIGARET

and passed to SUBSTRAT via the file CIGOUTl, to create the incident flux distribution.

TMPSUB2 created a prescribed flux distribution used by the program, using a simple

formula whose coeflSicients were obtained from its input file (magnitude of peak flux, a^,

Oy, prescribed velocity v). See equation (G2), Appendix G.

• SUBSTRAT now produces the file CIGINl, to be used by CIGARET.

• What was previously called "initial x-position of (the flux) peak" (on the substrate) is now the X = 0 end of the cigarette.

3. Documentation

Besides the present report, the internal documentation in CIGARET has been somewhat expanded over what was in CIG25, so as to make it easier to follow the program. Also, the NOMENCLATURE section

includes many of the symbols used in the program, which again facilitates understanding it, for anyone

who wishes to do so.

74
G. RESULTS

1. Sensitivity of Calculation

The temperature dependence of the reaction rate given by equations (108) and (1 10) is very strong. Thus, assume the peaked spatial dependence

T(x) « + Arexp[-(x-xJ^/o2] (132)

for T, where AT s T^^^ - T^. (See equations (C9) and (F2)). Then it is not difficult to show that the

width of the resulting reaction rate distribution is approximately a\8, where

T T

6 = ^ , (133)

Tm = Tniax. and Ta = E^/R is the activation temperature. T^ is of the order T^ ~ 22,000 K, while

Tnj w 1000 K. Hence 6 < < 1, and the reaction rate distribution R(x) is quite narrow. Hence the

gradient is high:

(134)

where R^, is the peak reaction rate. Integration of the reaction over a volume shows that we must also

have

R ~ (135)

Hence

|^«(ay6)-', (136)

dx

and the gradient thus rises rapidly with T^. Thus the grid size must be adequately small in order for the

numerics to adequately cope with this. In fact, with a grid size Ax = 0.2 mm, the program could not converge with T^s^ = 22,660 K, the measured value. The "brute force" approach would be to halve the

grid size to 100 microns. Not only would that quadruple the computational time to impractical levels,

it could not even be done with DOS-based PC's, because of inadequate memory space. It was therefore

necessary to use a lower, "model" value for T^^. To compensate, we also lowered the pre-exponential

factor such that the rate constant was approximatly correct in the temperature range of interest. The

program converges with T;^ = 15,000 K.

2. Velocity of Smolder Wave

Three 30-second runs were made with the input data shown in the sample input data. Section III.H.4, to obtain 90 seconds of free smolder by this hypothetical cigarette. Some of the results are shown in Figure

22. The upper and lower curves in that figure show the positions of the intersections of the 600 K

isotherms with the surface, as a function of time. Note that these isotherms enclose the reaction zone. Since the distance between them is still growing at t = 90 s, it is clear that a steady state has not yet been

reached.

75
0 10 20 30 40 50 60 70 80 90 100

t (sec)

Figure 22. Result of a sample run with CIGARET. The upper and lower

curves show the positions of the intersections of the 600 K

isotherms with the surface, as a function of time. The central

line is the locus of the midpoints between the curves.

76
Assuming that the isotherms are ^proximately symmetrical fore-and-aft, we have taken the midpoints

of the intersections of the 600 K isotherms at the surface to represent x^, "the" position of the smolder wave front. The fore-and-aft symmetry is not, in fact, perfect: the same procedure used for several

different isotherms at t = 60 s yields the values shown in Table 6.

Table 6. Position of Cigarette Smolder Front at 60 Seconds

as a Function of the Choice of Isotherm

Isotherm (K) Xo (mm)

600 6.4

700 6.6

800 6.7

Thus, as can be seen, although the "centroid" positions depend on which isotherm is chosen, the

differences are quite small, so that the 600 K isotherm is representative of the front.

The central curve shown in the figure is the mean value of the isotherms, which we may refer to as <x>goo- The points have been joined by a straight line, in fact, indicating that the smolder wave

(regression) velocity calculated this way is remarkably constant in the period shown, even though a steady

state has not (yet) been achieved. The slope of this line is V^oo = 5.36 nmi/min, a reasonable velocity.

The peak temperature in this calculation is not always on the axis. Although it is not plotted, the position

of the peak temperature is not a good indicator of the regression rate during these first 90 seconds. At

times it does not make any forward progress at all, for example. As the smolder continues, steady state would presumably be approached; when that eventually happens, of course the peak temperature will then

move at the same velocity.

3. Temperature Profiles

The temperature profiles at t = 60 s, from the run described above, are shown in Figure 23 as a function

of X, the distance from the cigarette tip. The upper curve is that on the axis; the lower one, the surface

temperature. Note that the axial temperature peak is at x = 3.2 mm, well forward of the surface peak,

at X = 5 mm. The interior peak evidently corresponds to the cone tip, and therefore it is quite

reasonable that the peak should be in front of the surface peak. Note, further, that the paper bum line

must be at 450 °C = 723 K, and that that occurs at x = 10.8 mm; hence the cone length is 10.8 - 3.2 = 7.6 mm, a not unreasonable length.

77
0 4 8 12 16 20 24 28

X (mm)

Figure 23. Longitudinal temperature distributions, along the axis and along

the surface, at t = 60 s, for a sample run with CIGARET.

78
H. USERS' GUTOE

1. Running CIGARET

CIGARET and two helper programs together with the SUBSTRAT program are available on a diskette

for IBM PC-compatible computers. The distribution diskette includes both executable and source code

for CIGARET and SUBSTRAT, executable code for the helper programs, and sample input files. The

CIGDATA program is used to prepaie data files for CIGARET. CIGDATA must be run on an IBM PCcompatible computer with VGA graphics. CIGARET requires a 486 class PC (no graphics needed) in

order to achieve satisfactory performance. On a 486/33 computer, 1 second of simulation time requires

about 6 minutes of clock time, a factor of over 300, for Ax = 0.2 mm (3172 nodes). Running the

program on a Silicon Graphics workstation, which is about 20 times faster than a 486/33 computer, is much more satisfactory.

The CIGARET source code (file CIGARET.FOR) can be compiled using any ANSI FORTRAN compiler.

All CIGARET input and output files are ASCII files. Therefore, CIGARET can be recompiled and run on a different computer, while still using CIGDATA and SUBSTRAT on a PC and transferring files between the computers. A diff'erent computer may allow CIGARET to execute faster.

The user should be sure to inspect the README file on die distribution diskette. One way to read this

file is to place the diskette in drive A: (or drive B:) and type

MORE <A:README or MORE <B:README

A permanent copy may be made with

PRINT <A:README

The README file contains a list of all files on the diskette, instructions for installing the necessary files on your hard disk, and information on any changes or additions to the program.

There should be at least 1 MB available on your hard disk. It is best to create a single subdirectory for

the executable programs and related data files. This will allow you to easily delete all the files related

to this program when you are finished with it. In general, keep allfiles in the current working directory.

Install the CIGARET program following the instructions in the README file. For example, from the

directory on the hard disk where you want CIGARET installed and with the diskette in A:, type

A:INSTALLA:

2. Input and Output Files

All input for CIGARET must be placed in a data file, whose contents are described in Section III.H.4.

This file is read as the "standard input stream," so CIGARET executes on MS-DOS and UNIX computers

by redirecting the input file. For example, begin a run by typing

79
CIGARET < CIGl.DAT

where "CIGl.DAT" is the name of your input data file. [There is a sample input file available on the

distribution diskette, TEST2.DAT]. As CIGARET completes each time step during the simulation, it displays the time (in seconds), the iteration number (K) within the current time step, and the time step

number (IT). This allows you to monitor the progress of the simulation. At first, it takes over 100

iterations per time step to converge; as time goes on, fewer and fewer iterations are required, until as few

as two or three sufllice.

Two files are written by CIGARET:

CIGOUTl provides cigarette surface temperature data for the SUBSTRAT program which are used to determine the incident heat flux on the substrate.

CIG0UT2 contains a dump of the cigarette field variables, which allows the program to be

restarted at any future time (see Section III.H.3).

3. Restarting a Run

CIGARET can be aborted in the usual way by pressing CTRL and C at the same time. The dump file CIG0UT2 will be produced.

A terminated run (either aborted or terminated in some other way) can be restarted by using the dump

file CIG0UT2. First, change its name to CIGIN2. Then, the input file (CIGl.DAT, in our example)

must be edited by changing the last data value in the file from 0 to 1; this alerts the program to use CIGIN2 to initialize the variables. Then type "CIGARET < CIGl.DAT" as before, and the execution

begins.

If a run is expected to be prohibitively long, it can be run in several steps. However, since a new copy

of CIGOUTl is produced every time CIGARET is run, it is necessary to save and then merge these files into a single CIGOUTl file. This is best done by renaming each file as it is created. For example, type

RENAME CIGOUTl CIGOUT.l,

RENAME CIGOUTl CIGOUT.2, etc.

Finally, combine the files by using the CMERGE program:

CMERGE CIGOUT.l CIGOUT.2 ... CIGOUT.N

The above command merges the specified files to create a new CIGOUTl.

80
4. Contents of the Data File

Line Variables Brief description

10

NR NZ RAD

PH PHI DP

Dl D2 D3

NR number of cells in radial direction

NZ number of cells in axial direction

RAD radius of cigarette [cm] (normalization; R)

PH void fraction

PHI total void fraction

DP diameter of pores [cm]

Dl tobacco density / gas density

D2 char density / tobacco density

AN ash density / tobacco density

SPHGS A FK SPHGS Cp(air) / Cp( tobacco)

D CN GK

B C E

TA TP EC

YA GAR GARP

ISP INP lEB

TSl TS2 ERR

A

FK

D

CN

GK

B

C

E

TA

TP

EC

YA

GAR

GARP

ISP

INP

lEB

TSl

TS2

ERR

1/ [density (solid) *Cp(so1)*r2 ] [cmK/J]

thermal conductivity of tobacco [W/cmK]

oxygen diffusion coefficient / R^ [s'^] mass of oxygen / mass of char consumed

thermal conductivity of gas [W/cmK]

pre-exponential factor * density [1/s]

heat of combustion / Cp( tobacco)

activation temperature [K]

ambient/ initial temperature (K)

paper ignition temperature [K]

cigarette surface emissivity, e^^

am±)ient oxygen mass fraction

mass transfer coefficient for air,/R:

combined coefficient (air+paper): R'*(Y,

Yg/R [1/9]

number of time steps

number of time steps between data outputs

number of time steps between energy balance checks

time step for first 50 steps [s]

time step for later iterations [s]

Gauss-Seidel convergence criterion

11 IDATA IDATA 1 = read CIGIN2 to restart CIGARET

81
Sample Data File

(Note: the number of spaces between inputs on a line is arbitrary.)

21 151 0.40

0.65 0.85 0.0575

627. 0.341 0.13

0.777 6.5 3.16E-3

0.7 1.6428 4.514E-4

1.1085E9 13461.5 1.50E4

293.15 723.15 0.73

0.232 8.225 2.725

10000 200 80000

0.005 0.005 0.001

0

5. Producing a Data File

The CIGARET input data file can be created with any ASCII line editor. CIGDATA creates these data

files interactively and thus uses certain commands which restrict its operation to IBM PC-compatible

computers. It includes some checking of the input data. CIGDATA is especially useful for creating a data file which is only slightly difl'erent from another data file. This is useful in performing the

parametric studies for which CIGARET was designed.

There are two special files in CIGARET to help the user with CIGDATA:

The help file, CIGDATA.HLP, contains the text of the interactive help messages. Help is activated by pressing the Fl function key. If the help file is not available in the current working

directory, no interactive help will be available.

The configuration file, CIGDATA. CFG, sets the colors of the display. The file included on the

distribution diskette assumes that a standard VGA monitor is being used. If the configuration file

is not in the current working directory, a set of default colors will be used. A new configuration

file can be made by using the MAKECFGT program. See the README file for instructions.

The operation of CIGDATA is explained on the following pages, which show the messages and input screens which will appear as the program is run. After reading through these pages, try using CIGDATA

with one of the sample data files. Begin the program by typing CIGDATA. Abort the program by

pressing CTRL and C simultaneously.

82
Data input screens from CIGDATA:

Screen 1: (primary menu)

CIGARET data preparation:

# File information

# Geometry data

# Tobacco data

# Pyrolysis data

# Gas data

# Boundary data

# Simulation control

# Exit data preparation

Use cursor keys to move between menu selections.

Press ENTER to activate the menu selection at the X.

Press ESC to return from a selection. Press Fl for help.

Screen 2:

File data:

Neune of old data file:

Name of new data file:

press ESC when done; press Fl for help

Screen 3:

Geometry data:

Radius of cigarette:

Number of cells in radial direction:

Number of cells in axial direction:

press ESC when done; press Fl for help

[mm]

83
Screen 4:

Tobacco data: press ESC when done; press Fl for help.

Density of tobacco:

Thermal conductivity:

Specific heat:

Void fraction:

Total void fraction:

Diameter of pores:

Screen 5:

Reaction data:

Pre-exponential factor:

Activation temperature:

Heat of reaction:

Oxygen mass / char mass:

Density of char:

Density of ash:

press ESC when done; press Fl for help.

[m^3/kg s]

[K]

[kJ/kg]

[kg/m-3]

[kg/m-3)

See the explanatory note on the next page re "Activation temperature."

Screen 6:

Gas data: press ESC when done; press Fl for help.

Thermal conductivity of the gas: [W/m K]

Specific heat of the gas: ^^HHSU^ [kJ/kg K]

Oxygen diffusion coefficient: [m^2/s]

84
Screen 7:

Boundary data: press ESC when done; press Fl for help.

Initial/ambient temperature: ^^^^J^^^^ [°C]

Paper ignition temperature: ^^^^Jjl)^^

Cigarette emissivity:

Ambient oxygen mass fraction:

Mass transfer coefficient through boundary layer:

Mass transfer coefficient through virgin paper:

[m/s]

[m/s]

Screen 8:

Simulation control: press ESC when done; press Fl for help.

Total number of time steps: ^^^^M

Steps between outputs to substrate:

Time step for first 50 steps:

Time step for later iterations:

Gauss-Seidel convergence criterion:

(S) begin simulation with default start-up conditions

{ft) read CIGIN2 to resume a simulation

Activation Temperature: As in TMPSUB2 (and therefore in SUBSTRAT) , we use the

activation temperature of the reaction, defined as:

(where E^ is the activation energy and R is the universal gas constant) ^ as the

input parameter, rather than E^ itself.

85
I. SUMMARY OF SECTION III

The dynamics of a freely smoldering cigarette have been discussed. Similarly, what happens when it is resting on a substrate has been considered. Some of the mathematical models designed to simulate a smoldering cigarette have been outlined. Then the equations that must be satisfied by a cylindrically

synmietric, homogeneous model of a cigarette quietly (freely) smoldering in air have been set down,

within some specified simplifying assumptions (Section III.D.l). The effects of the paper wrapping,

which pyrolyses away, are included.

The present program, CIGARET, is the outgrowth of an earlier model, CIG25, which employed some

simplifications to the equation; principally, the neglect of pyrolysis (which has relatively little consequence in terms of the energy balance). CIGARET is a much-improved program in a number of

ways which are described in detail in Section III.F. The numerical method used for solving the equations

is described, and some of the problems discussed, in Section III.E. A primer on the use of CIGARET

is given in Section III.H.

Some of the results obtained from using CIGARET are given in Section ni.G. One remarkable result

is that after ignition, the velocity of the resulting smolder wave quickly becomes constant, long before

a steady state is achieved. The region in which char oxidation is proceeding rapidly (the glowing region)

is correctly calculated to be conoid-shaped, and of the correct length.

86
IV. SIMULATING A BURNING CIGARETTE ON AN IGNITABLE

SUBSTRATE

A. Introduction

The process which we are simulating is the heating and possible ignition of a substrate by a cigarette. In order to do that, the models SUBSTRAT and CIGARET developed in Sections II and III must be used

in conjuction with each other. In this Section, the effects which take place when the lit cigarette and the

substrate are in intimate contact are discussed. Next, it is shown in Section IV. C how the two programs

are to be used in tandem in order to simulate the interaction. Finally, the effects are then calculated in

detail. It is not essential that the reader read these last subsections in order to be able to use the program

intelligently.

B. Qualitative Description

Consider the effects of a cigarette on a horizontal substrate that might affect the ignition of the substrate:

• The principal effect is the heating of the substrate by the hot cigarette coal.

• If any oxidative reaction takes place in the substrate, the cigarette competes with the

substrate for oxygen.

• Some of the water vapor and tar emitted by the cigarette may recondense on the surface

of the substrate and change its thermal characteristics.

• The cigarette affects the boundary conditions on the substrate: it interferes with

convective cooling over the entire length of the cigarette.

Conversely, the horizontal substrate will influence how the cigarette smolders and thereby affect the

likelihood that the cigarette will ignite it:

• The cold substrate initially provides a substantial conductive heat sink to the cigarette,

• If the substrate eventually undergoes exothermic reactions, it will heat the cigarette,

instead.

• The substrate obstructs access of oxygen to the cigarette.

• If the substrate begins to react with oxygen, then the oxygen depletion for the cigarette

would become still more severe.

These effects fall into two categories: gas transport effects and thermal effects. Consider the former first. There is a boundary layer surrounding the cigarette, within which the oxygen concentration drops as the

surface is approached. Oxygen is supplied to the combustion site by the ambient via diffusion, as indicated by the arrows in Figure 23. The proximity of the substrate reduces the availability of oxygen

to the cigarette from below. This will tend to inhibit the smoldering rate. This effect can be simulated

by assuming that y^ < 0.232; CIGARET is given the information through the input file. See Section

IV.D.4.

87
Thermally, the colder substrate cools the cigarette near the line of contact with the substrate, further

inhibiting smoldering. This sink is not constant: as time progresses, the substrate heats up, and its cooling effect becomes progressively weaker. Radiative losses from the cigarette are decreased both

because the cigarette is cooler and because one of the "targets" for this heat is warming and re-radiating

back to the cigarette. If the heat loss from the cigarette is sufficiently high, the cigarette goes out. This

effect is schematically indicated as curve A in Figure 21. If the substrate is not inert and exothermic

reactions take place, cooling may be replaced by heating; this is sketched as curve B in the figure. This

effect is also given in the input file for the cigarette, through the surface boundary conditions, as discussed in Section IV.D.

When both the cigarette and the substrate are reacting with oxygen, they are concurrently (a) generating

heat, which accelerates their heat production, and (b) competing for oxygen, which retards their heat

production. The modeling of effect (a) has been extensively described in Sections II and III; the modeling

of (b) is not treated in this study.

C. Use of the Two Programs

It is important to understand how CIGARET and SUBSTRAT are used to obtain the interaction between

the smoldering cigarette and the substrate with which it is in contact. The cigarette-substrate interaction

is obtained through the boundary conditions which apply to each of the two models.

A fiiU simulation requires an iterative use of CIGARET with SUBSTRAT (Section IV). Before starting, one should be sure there is no file named CIGINl. First, run CIGARET (see Section III.H), generating

CIGOUTl. The latter is now read by SUBSTRAT by typing

SUBSTRAT < INPUTFILE

(see Appendix B). [This was done before for TMPSUB2, as well.]

SUBSTRAT now produces a new input file, CIGINl, for the CIGARET program. This file contains

substrate surface temperature data which become part of the boundary conditions for the cigarette. Then

run CIGARET, which will read CIGINl and create a new CIGOUTl. This simulation now includes an

estimate of the influence of the substrate on the cigarette. Then run SUBSTRAT again, and so forth until convergence is obtained. The input datafiles produced by the userfor CIGARET andfor SUBSTRAT (as

distinctfrom thefiles produced by the programs themselves) must remain the same during these iterations.

Since the principal effect of the cigarette is to heat the substrate, how to find the heat flux from the

cigarette to the substrate will now be considered in detail. Then, the reciprocal effects of the substrate on the cigarette are discussed.

88
D. DETAILED CALCULATIONS OF INTERACTIONS

1. Conduction Flux to Substrate

We begin by calculating the heating flux delivered to the substrate by the cigarette. The substrate starts out at the ambient temperature and is in fair thermal contact with the cigarette. As shown in Figure 24,

the measured flux is the initial (and momentary) flux from the cigarette to the cold substrate at the very

small area of contact between them, along the LC Oine of contact). The peak measured flux shown is about 5.6 W/cm^. Since the measured coal surface temperature there was about 550 °C, the (net)

blackbody radiation to the ambient is 2.56 W/cm^. Muramatsu measured = 0.73; hence the actual radiation loss (assuming a grey body) is (^^ad =1-9 W/cm^. This leaves 5.6 - 1.9 = 3.7 W/cm^ loss

rate via conduction. With a surface temperature of 550 °C and an ambient temperature of 20 or 25 °C, one infers that the eff^ective heat transfer coefficient is 71 W/m^K, about seven times what it is in air. This large flux (5.6 W/cm^) heats the substrate rapidly, so that the net flux, given by the first part of

equation (137), falls rapidly to lower values.

The flux distribution <^(x) shown in Figure 24 is that along a narrow region (1.5 mm wide) about the line

of contact; that is, it is 0(x,O). In order to find the heating of the substrate, we must also know the transverse and time dependences of the flux: <t> = 0(x,y,t).

It has been assumed that the initial substrate surface temperature is the ambient temperature, i.e. , Tg(x,0) = T^. The heating flux from the cigarette to the substrate along the CL, <l>^(\,0,t), is given by experiment, as in Figure 24; the flux from the (resulting) hot substrate, 0s(x,t), is calculated by the substrate

program.

Moving away from the line of contact, the temperatures rapidly approach ambient. Therefore, the

convective heat transfer eff'ects of the cigarette diminish and the heat transfer coefficient asymptotically

falls to its ambient value. A reasonable approximation to how the net convective flux to the substrate

must vary with y is

where h^. is the heat transfer coefficient for convective heat loss from the cigarette to the substrate at the

point at the line of contact corresponding to the peak flux (about 71 W/m^K). h^, is the heat transfer

coefficient for convective heat loss from the surface to the ambient (about 10 W/m^K), and Q^Cy) is a

fraction indicating the convective inffuence of the cigarette at the distance y. Note that the net convective

flux "into" the surface can become negative. If the cigarette extinguishes, for example, the fluxes from

the cigarette disappear, and T^. - Tg is asymptotically replaced by T^ - Tg.

We can break this flux into two parts, one of which is independent of the substrate temperature, and the

other of the cigarette temperature. Thus, the cigarette "emits" the convective flux

^ ho\^T^ - T,(x,y)][l - Q^(y)]

(137)

(138)

(139)

89
Figure 24. Flux emitted by a cigarette toward the substrate,

along the contact line. The dashed curves

indicate the probable errors (variance).

90
It is shown in Appendix G that QJy) can be approximated by a Gaussian. It follows that the effective

heat transfer coefficient in equation (139) becomes

where

ho = background « 10 W/m^K and Ah = h^ - h^ « 61 W/m^K.

What is to be used for ffy in equation (140) is discussed in Section rV.D.2, just below.

2. Radiative Flux to the Substrate

The radiation exchange experienced by the substrate at point (x,y) is

K..r = e,Qa(r* - r/) - (1 -Q)€,0(t/ - r/), (Hi)

where the x,y arguments have again been left out for clarity, and where

(142)

' 1 -a

is the effective emission coefficient. The radiation flux is treated in the same way as the convective flux;

tiiat is, this is split into an effective radiation flux from cigarette to substrate,

*.„, = .,Qa(r/-r;), (143)

while the substrate radiates away at the rate W= [^^" *€,(i-Q)]a(r/-r;) (144)

In Appendix 5-D of Gann et al. (1988), it was shown that the variance of the Gaussian approximation

to the convective distribution is about 0.57R. In Appendix H of the current work, it is found that fi is reasonably well-approximated by a Gaussian also, of variance s. In Appendix G, it is shown that the joint

convective-plus-radiative flux distribution is also approximated by a single Gaussian. We therefore use

the approximation

Qc(y) "(y) = exp(-yVa5) (I'^S)

(See equation (H13)). Since we found above that Oy = 3.2 mm for R = 4nmi, we may assume, in

general, that the net variance is 0.8R:

Oy = 0.8/? (146)

(Also see equation (H14)).

91
3. Complete Flux

The cigarette, substrate, and net fluxes are the sums of the convective and radiative contributions.

Adding together the terms from the two previous subsections, we can, for the general case, write this in the form

Here the flux at (x,y) from the cigarette is

= K^^iT^ - - €,Q a(r/ - r/) (i48)

while the loss from the substrate to the ambient surroundings at (x,y) is

Kt(^^y) = K(T, - M€,Q - €,(1 - Q)] o(r; - rj*) (m

where the x-dependencies of T,. and are not shown explicitly in the terms on the right, for the sake

of brevity. These model fluxes hold for 0 < x ^ L. <f)^^^ is the model flux from the cigarette; <})^^^^ is the model flux loss from the substrate. We note that this separation does not quite succeed in giving a flux <j>^^^ which depends only on

cigarette/air properties and another which depends only on substrate/air properties. However, (i>^-^^ does

not depend (as it must not) on Ts(x,y,t), and <j)^^^^^ does not depend on Tc(x,y,t). The separation into 0^.jg

and </)sub is thus good enough that one can run the CIGARET and SUBSTRAT programs independently.

One obtains the efl'ects of the cigarette on the substrate from using, as input to SUBSTRAT, the flux <l>^^^ which is independent of Tg. On the other hand, by using 0^^^' one gets the surface losses of the substrate

without having to know T^(x). Of course the substrate "knows" about the cigarette through the latter's flux, <^cig» which is entered as an input to CIGARET. The substrate temperature equally influences the

cigarette, but by a slightly difi'erent mode; that is discussed in the next Section.

4. Effects of the Substrate on the Cigarette

We next consider the inverse problem. The energy losses of the cigarette will be afl'ected by the presence

of the surface. Those efl'ects are estimated in this Section. Moreover, since CIGARET is a model of

afreely smoldering cigarette, the only way we can determine the eff'ect of the substrate on it is to express

the cigarette loss rate so that it appears to be simply losing energy to the air. This can be done by using an efl'ective heat transfer coefficient and an eff'ective surface emissivity, such that the presence of the

substrate is correctly, if "covertly," taken into account. This is in fact possible to do, as is shown below.

Thermal Eff'ects: Convection. It is shown in Appendix F that the convective loss of the cigarette can be

expressed as:

Kcon « h^T^-TJ, (150)

92
where

h* - 21.22 - 10.47 - 9.02

c c.min T - T

W/m^K (FIO)

In deriving this, it was assumed that there is no longitudinal dependence of temperature for either the

cigarette or the surface. We now relax this assumption, and recognize that there are longitudinal

variations for both T^, and T^. To be consistent with the approach outlined below equation (Hll), we

assume equation (FIG) to be valid at every x, independent of adjacent regions. That immediately gives

the x-dependence of h* to some approximation. This will likely be a weak dependence. Consider the

middle term in equation (FIG). Along a constant-y line (e.g., y = y^), the smaller T^(x) is, the smaller

will be Tg(x,yQ), since that substrate surface temperature was produced by T^.. That suggests that the

ratio (Tg - TJ/(Tg - TJ in equation (FIG) will not vary strongly with x.

In order to use equation (FIG), we must know Tc,m\n- '^'s temperature is defined in Appendix F; it is the lowest temperature around the cigarette circumference, and lies at the contact line. This varies with

time, as well as with x; hence h* = h*(x,t). Also t, as defined by equation (F8), is a function of t. The resulting T^, n,in(t) is found in Appendix I. The result is that the temperature difference between the

top and the bottom of the cigarette is 6T^^, given by

6 = 0.62 r f (15) max where _ T + T

f = " ' (II)

1 +r

and where r is the ratio of the thermal inertias of the materials:

r = (12)

We assume that this minimum temperature is reached in two seconds. It is indicated in Appendices F

and I that this is the drop at the flux peak. The drop will be smaller away from the peak, and we may

expect something like

T,(x^,t)-T^

(151)

to hold.

It has been observed that it takes quite a long time for the cigarette to "recuperate" from this drop in (lower) surface temperature, assuming it does not, in fact, extinguish. Again for the sake of simplicity,

we will take the recovery period to be 2 minutes.

Suppose the cigarette is dropped onto the substrate at time t^. Then collecting all these observations, we

obtain the estimate

93
( 122 +t-t\

120

t+2<t<t+122 ^^^^^

o o

0 t > t+122 O where

T(x,0 = T^{Q = nl2,x,t) - r^(e=0,x,0 (F8)

With these specifications, equation (FIO) can be implemented.

Thermal Effects: Radiation . We next determine how the radiation losses of the cigarette are modified

by the presence of the substrate. The radiation loss for the free cigarette at any point on its surface is given by equation (143), with replaced by e^. For the cigarette on the substrate, along the contact line,

it is given by equation (H16). Notwithstanding the arguments made earlier about the difficulties in calculating the relevant radiation view factors, this was done, making some substantial simplifications.

We assume that the radiation flux emanating from the substrate is

^.(y) - 4), ^ A<|)exp(-y2/p2) (153)

where

A(J) = €^oT,\y=0) - 4), (154)

and

. (155)

The x-dependence will be discussed in a moment. A certain fraction of that flux reaches the cigarette.

Upon finding that fraction and integrating over y, it is found that the radiation loss of the cigarette, when

in the presence of a substrate, is cut down from the usual expression to

^c.r ' ^A^T,' - ^A<1)] - €,4), (156)

where T^, is assumed to be independent of 6. (See Appendix F). The approximate factor 0.7 is the result of a large number of simplifications and numerical integrations. One of the assumptions is that p = (Ty.

In the T3(0) = T^ limit, equation (156) is approximately correct, since aT^"^ > > 0^, the deviation from

exactitude is negligible. For higher temperatures T^(0), the flux <j)^^ decreases, as it should, and is still a substantial fraction of crT^,'* for the upper limit, T3(0) = T^..

Equation (156) can be written as Jar/-., (157)

where

94
1 - 0.7A<t>

(158)

Finally, this must be generalized to include the dependence on x and t. Just as was done for the

convective part, the simplest (although crude) way to do this is to generalize e^* by fiat to ej.*(x,t), where

eUx,t) = e ^ 0.7A4)(j:,0

e,ar/(x,f)

(159)

and

A4)(x.O = e^aT*(x,0,t) - <{), (160)

CIGARET can then be run without explicitly introducing the substrate temperature by using h*(x,t) and

ec*(x,t) as the model heat transfer coefficient and cigarette emissivity. These are easily computed at each

time step, in CIGARET.

In summary, explicit expressions for the transverse dependence of the cigarette and substrate fluxes have

been found, which supplement the longitudinal dependencies. Thus, if the longitudinal dependence on

the cigarette surface is known, the flux <^(x,y,t) can be found.

Fluid Flow Efl'ects . Quite independently of any possible "competition" from the substrate for oxygen

(that is, for oxidative reactions), access to air is somewhat restricted from below because of the presence

of the substrate. Thus, the efl'ect of having the substrate there is equivalent to limiting the indrawnoxygen rate to what it would be in the open, but with y^ < 0.23. This eff'ect is incorporated in running

CIGARET through the appropriate value for y^^ being supplied by the user in the input file.

In order to see what a small reduction in y^ produces, CIGARET was run with the same input parameters

as yielded the results shown in Figures 22 and 23, but one: y^ was taken as 0.21. The result, however, was that the simulated cigarette extinguished after a few seconds. [The criterion used to declare the

cigarette "extinguished" is that the peak surface temperature falls below 700 K (427 °C)]. Yet, when

the cigarette surface temperature history of the y^ = 0.23 run was translated into a flux, and the substrate

described in Section II. E was exposed to this flux, the substrate ignited in 22 seconds. Clearly, these are contradictory results, and the efl'ective value of y^ to be used for this cigarette lying on this substrate

should be between 0.21 and 0.232. It will probably not be very diff'erent for other cigarette/substrate

combinations.

V. ACKNOWLEDGEMENTS

The authors wish to thank Dr. Thomas J. Ohlemiller for the many insightful discussions on various

aspects of the work, Dr. Richard G. Gann for his careful and helpftil reading of the manuscript, and the

Consumer Products Safety Commission for sponsoring this work.

95
VI. REFERENCES

Ames, W.F. (1969) Numerical Methodsfor Partial Differential Equations, Barnes and Noble, New York.

Baker, R. (1975) "Temperature Variation Within a Cigarette Combustion Coal During the Smoking

Cycle," High Temperature Science 7, p. 236.

Baker, R.R., and Crellin, R.A. (1977) "The Diffusion of Carbon Monoxide Out of Cigarettes," Beitrage zur Tabakforschung 9, No. 3 (October), p. 131.

Belytschko, T. (1983) "An Overview of Semidiscretization and Time Integration Procedures," Chapter

1 in Computational Methodsfor Transient Analysis, edited by T. Belytschko and T.J.R. Hughes, NorthHolland.

Bird, R.B., Stewart, W.E., and Lightfoot, E.N. (1960) Transport Phenomena, Wiley, NY.

Brandrup, J., and Immergut, E.H., eds. (1990) Polymer Handbook, 3rd ed., Wiley, NY.

Carslaw, H.S., and Jaeger, J.C. (1959) Conduaion ofHeat in Solids, 2nd ed; Oxford University Press.

Childs, G.E., Ericks, L.J., and Powell, R.L. (1973) "Thermal Conductivity of Solids at Room

Temperature and Below," NBS Monograph 131, National Bureau of Standards, Gaithersburg, MD.

Clausing, A.M. (1969) "Numerical Methods in Heat Transfer" in Advanced Heat Transfer, Ed. B.T,

Chao, University of Illinois Press, Urbana, IL.

Cohen, N. (1961) "Boundary-Layer Similar Solutions and Correlation Equations for a Laminar Heat

Transfer Distribution in Equilibrium Air and Velocities up to 41,000 Feet per Second," NASA Technical

Report R-118.

Conte, S.D. and de Boor, C. (1972) Elementary Numerical Analysis, McGraw-Hill, NY.

Croft, D.R., and Lilley, D.G. (1977) Heat Transfer Calculations Using Finite Difference Equations,

Applied Science Publishers Ltd., London.

Egerton, A., Gugan, K., and Weinberg, F.J. (1963) "The Mechanism of Smoldering in Cigarettes,"

Combustion and Flame 1, p.63.

Gann, R.G., Harris, R.H. Jr, Krasny, J.F., Levine, R.S., Mitler, H.E., and Ohlemiller, T.J. (1988)

"The Effect of Cigarette Characteristics on the Ignition of Soft Furnishings," NBS Technical Note 1241,

National Bureau of Standards, Gaithersburg, MD.

Gugan, K. (1966) "Natural Smoulder in Cigarettes," Combustion and Flame 10, p. 161.

Hilsenrath, J., Beckett, C.W., Benedict, W.S., Fano, L., Hoge, H.J., Masi, J.F., Nuttall, R.L.,

Touloukian, Y.S., and Woolley, H.W. (1955) Tables of Thermal Properties of Gases, NBS Circular 564,

National Bureau of Standards, Gaithersburg, MD.

Holman, J.P. (1981) Heat Transfer (5th ed.), McGraw-Hill, New York.

96
Kaka?, S., Shah, R., and Aung, W., eds. (1987) Handbook of Single-Phase Convective Heat TYansfer,

Wiley/Interscience, New York.

Kashiwagi, T., and Nambu, H. (1992) "Global Kinetic Constants for Thermal Oxidative Degradation of

a Cellulosic Paper," Combustion and Flame 88 p.345.

Kunii, D. (1961) Kagaku Kogaku [Chem. Eng. (Japan)] 25 p.891.

Larkln, B.K. (1964) "Some Stable Explicit Difference Approximations to the Diffusion Equation,"

Mathematics of Computation 18 p. 196.

Lawson, J.R. (1993) National Institute of Standards and Technology, unpublished data.

Lendvay, A.T., and Laszlo, T.S. (1974) "Cigarette Peak Coal Temperature Measurements," Beitrdgezur

TabaJtforschung 7, p. 63.

Miller, A.L. (1991) "Where There's Smoking There's Fire," NFPA Journal, p.86.

Mitler, H.E., and Davis, W.D. (1987) "Computer Model of a Smoldering Cigarette," Annual Conference on Fire Research, National Bureau of Standards, Gaithersburg, MD, Nov. 1987.

Mitler, H.E., and Walton, G. (1992) "A Computer Model of the Smoldering Ignition of Furniture,"

NISTIR 4973, National Institute of Standards and Technology, Gaithersburg, MD.

Peyret, R., and Taylor, T.D. (1983), Computational Methods for Fluid Flow, Springer-Verlag, New

York.

Moussa, N.A., Toong, T.Y., and Garris, C.A. (1977) "Mechanism of Smoldering of Cellulosic

Materials," Sixteenth Symposium (International) on Combustion, The Combustion Institute, Pittsburgh,

PA, p. 1447.

Muramatsu, M., Umemura, S., and Okada, T. (1979) "A Mathematical Model of Evaporation-Pyrolysis

Processes Inside a Naturally Smoldering Cigarette, Combustion and Flame 36, p. 245.

Muramatsu, M. (1981) "Study of Transport Phenomena Which Occur During Unforced Smoldering of

Cigarettes," Res. Report No. 123 of the Japan Tobacco and Salt Monopoly (in Japanese).

Ohlemiller, T.J. (1985) "Modeling of Smoldering Combustion Propagation," Progress in Energy and

Combustion Science 11, p. 277.

Ohlemiller, T.J. (1991) "Smoldering Combustion Propagation on Solid Wood," in Fire Safety Science,

Proceedings of the Third International Symposium (Eds., G. Cox and B. Langford); Elsevier Publishing

Co., London, p. 565.

Ohlemiller, T.J. (1991) unpublished results.

97
Ohlemiller, T.J., Villa, K., Braun, E., Eberhardt, K.R., Harris Jr., R.H., Lawson, J.R., and Gann, R.G,

(1993) "Test Methods for Quantifying the Propensity of Cigarettes to Ignite Soft Furnishings," Special

Publication 851, National Institute of Standards and Technology, Gaithersburg, MD.

Ozisik, M.N. (1980) Heat Conduction, Wiley-Interscience, New York, 1980.

Parker, W.J. (1985) "Prediction of the Heat Release Rate of Wood," Fire Safety Science - Proceedings

of the First International Symposium (Eds., C.E. Grant and P.J. Pagni); Hemisphere Publishing Corp,

p.207.

Parker, W.J. (1988) "Prediction of the Heat Release Rate of Wood," Ph.D. Thesis, George Washington

University, Washington, DC.

Quintiere, J.G., Harkleroad, M.F., and Walton, W.D. (1983) "Measurement of Material Flame Spread

Properties," Combustion Science and Technology 32 p. 67,

Quintiere, J.G., and Harkleroad, M.F. (1985) "New Concepts for Measuring Flame Spread Properties,"

Fire Safety Science and Engineering, ASTM STP 882, T.Z. Harmathy, ed., ASTM, Philadelphia, PA,

p.239.

Quintiere, J.G. (1988) "The Application of Flame Spread Theory to Predict Material Performance," J. of Res. ofNat'lBur. of Standards 93{l), p.61.

Raznjevic, K. (1976) Handbook of Thermodynamic Tables and Charts, Hemisphere Pub. Co., New York.

Samfield, M. (1986) private communication.

Sandusky, H.W. (1976) "A Computer-Simulated Cigarette Model for Use in the Development of Less

Hazardous Cigarettes," Ph.D. Thesis, Princeton University, Princeton, NJ.

Schneider, P.J. (1973) "Conduction," Handbook of Heat Transfer, 1st Ed., Eds. W.M. Rohsenow, &

J. P. Hartnett, McGraw-Hill, New York.

Siegel, R., and Howell, J.R. (1981) Thermal Radiation Heat Transfer, 2nd ed., McGraw-Hill, NY.

Sparrow, E.M., and Cess, R.D. (1978) Radiation Heat Transfer, McGraw-Hill, NY.

Summerfield, M., Ohlemiller, T.J., and Sandusky, H.W. (1978) "A Thermophysical Mathematical Model

of Steady-Draw Smoking and Predictions of Overall Cigarette Behavior," Combustion and Flame 33,

p.263.

Szekely, J., Evans, J., and Sohn, H. (1976) Gas-Solid Reactions, Academic Press, New York.

Thomas, P.H. (1957) "Some Conduction Problems in the Heating of Small Areas on Large Solids,"

Quart Journ. Mech. and Applied Math., Vol. X, Pt. 4.

Torrance, K.E. (1985) "Numerical Methods in Heat Transfer," Handbook of Heat Transfer Fundamentals, 2nd Ed., Eds. W.M. Rohsenow, J. P. Hartnett, and E.N. Gacic, McGraw-Hill, New York.

98
Touloukian, Y.S., Powell, R.W., Ho, C.Y., and Klemens, P.G. (1970) Thermophysical Properties of

Matter, Vol.2, Thermal Conductivity, Nonmetallic Solids, IFI/Plenum.

Treybal, R.E. (1955) Mass-Transfer Operations, McGraw-Hill, New York.

Weast, R.C., Ed. (1976) Handbook of Chemistry and Physics, CRC Press, Cleveland, OH.

Wichmann, I.S. (1991) "On the Use of Operator-Splitting Methods for the Equations of Combustion,"

Combustion and Flame 83, p. 240.

99

APPENDIX A

CONDUCTION ALGORITHM TESTS

A computer program called TEMPSUB was developed as part of the earlier investigation into the ignition

of furnishings by smoldering cigarettes (Garm et al., 1988). This program modeled heat transfer in

furniture, or in a "substrate," using a simple finite difference approximation (FDA) for a homogeneous

substrate with uniform and constant properties. The research indicated that this program would have to be expanded to include a two-layer model (fabric + padding), pyrolysis of each layer, an asymmetric flux input, and a variable grid.

These features have been implemented in SUBSTRAT by using a slightly different approach to the FDA

than was used in TEMPSUB. The original approach was to convert the differential equation for heat

transfer into an FDA by a Taylor's series approximation. The new approach is based on tiie conservation

of energy within a control volume. It is based on physical reasoning and is usually easy to apply. It is most useful for variable grids, convective boundary conditions, odd-shaped regions, etc. It is more

difficult to obtain accuracy estimates for the control volume approach than for the Taylor's series

approach. For simple cases involving uniform grids and homogeneous materials the two approaches lead

to identical FDA's. See Torrance (1985) or Croft and Lilley (1977) for further details.

Since there is a considerable increase in the desired capabilities of the upgraded program, it was decided

to develop a program which would allow extensive testing of the FDA. This program is called CTEST3

(Conduction TEST - 3 dimensional) which has the ability to model simple boundary conditions (constant

temperature, heat flux, or convection coefficient) on any of the six faces of the region.

The FDA used in CTEST3, the explicit Euler method, has been checked against several heat transfer

problems which have analytic solutions. The first few tests involve various combinations of constant temperature, heat flux, and convection coefficient boundary conditions with analytic solutions from the

classic text by Carslaw and Jaeger (1959).

Several of these analytic solutions involve the error function which is defined by

X (Al)

so that erf(0) = 0,

erf(oo) = 1,

and erf(-x) = - erf(x).

The complementary error function is also used. It is defined as

(A2)

so that erfc(0) = 1,

and erfc(oo) = 0.

A-1
Repeated integrals of the error function are also useful in conduction problems. These are defined by

the recursive relationships

i'^erfcW = I

i''-^erfc(x)dx n=l,2,... (A3)

X

or

i<*erfc(x) = &rfc(x) (A3a)

ierfcW = i^erfc(x) = ;cerfc(x) (A3b)

Ini'-erfcCx) = i^-^erfcW - 2jci''-Wc(x) n=2,3,... (A3c)

See Appendix II of Carslaw and Jaeger (1959) for further details on error functions. Computer

subroutines were written implementing these functions for the computation of the analytic solutions of

the heat transfer tests.

Test 1: One-Dimensional Steady-State Conduction

CTEST3 has been tested for steady-state conduction with constant thermal properties and uniform grid

spacing. This test sets opposite faces on a cubic region to different temperatures and makes the remaining

faces adiabatic. After a sufficient number of time steps the temperature within the region should vary

linearly from the hot to the cold face:

Tlx) = TXO) + -[IlL) - IXO)] (A4)

This has been confirmed in all three directions.

Test 2: One-Dimensional Transient Conduction, Constant Heat Flux Boundary Condition

This test was used in the development of the original substrate model. The analytic case for a constant

flux involves a homogeneous solid occupying the semi-infinite region x > 0. The solid is initially at zero temperature throughout. At time t = 0, a constant heat flux, q, is applied to the x = 0 surface. The

temperature within the region is given by Carslaw and Jaeger (1959), p. 75, equation (6):

7Xx,0 = i^^erfc

K

X

(A5)

Preliminary testing (again using a uniform grid and constant thermal properties) indicated that as At and

Ax decreased, there was a uniform approach to the analytic solution. As for the accuracy to be expected,

for Ax = 1 mm and At = 0.5 s, the error (after the first three time steps) was < 1 percent. We note

that the results of test did not agree with results from the original TEMPSUB model. Further

investigation indicated an error in the TEMPSUB boundary conditions subroutine. Correcting this error brought results from the two programs into complete agreement.

A-2
Tests show that reducing the grid size along with corresponding reduction of the time step cause the FDA

solution to approach the analytic solution. Therefore, the FDA is consistent. Reducing the time step

without changing the grid size does not improve the accuracy of the solution. In fact, it is best to operate

as close to the stability limit as possible for botii accuracy and execution time. Use of the variable grid

gives results consistent with the uniform grid at the surface. The error in the calculated surface

temperature goes down as time increases.

Test 3: One-Dimensional Transient Conduction, Constant Convection Coefficient Boundary

Condition

This test represents a slab of a homogeneous solid of diickness 2L in the x direction and infinite extent

in the y and z directions which is initially at unit temperature throughout. At time t = 0 the temperature

of the fluid on both sides of the slab is changed to zero and heat is convected from the slab through a

constant convection coefficient. Because of symmetry this problem is equivalent to a slab of thickness

L with one adiabatic surface at x = 0 and a convective surface at x = L. The temperature within the

region is given by Carslaw and Jaeger (1959), p. 122, equation (12):

" 2Bcos(6„x/I)

T(.x,t) = T "

where B = hL/<c (Biot number),

F = at/L? (Fourier number), and

6n are the solutions of the transcendental equation

(A6)

6„tan6„ = B.

In order to check the calculation of this complicated analytic solution, the solution to a related problem was also computed. This is the temperature in a semi-infinite slab with the convective boundary condition

(Carslaw and Jaeger (1959), p. 72, equation (5)):

T(x,t) = erfc - exp

(h

—(x+aht/y.)

K

erfc (A7)

where x is now the distance from the convective surface into the region. There is good agreement

between the FDA and analytic solutions. Again, the error in the calculated temperature goes down as time increases.

Test 4: Three-Dimensional Transient Conduction, Constant Surface Temperature Boundary

Condition

This test represents a block of a homogeneous solid in the region defined by -a < x < a, -b < y < b,

and -c < z < c which is initially at unit temperature throughout. At time t = 0 the temperatures of the

surfaces of the block are reduced to zero, and the block begins to cool. The temperature within the

region is given by Carslaw and Jaeger (1959), p. 184, equation (5):

A-3
64 (-1)i+JII+ll Uti^Ti (2/+l)(2m+l)(2n+l)

cos

(2/+l)7IJC

2a

cos (2m+l)Tty

2^> cos (2n+l)Kz

2c (A8)

(2/+l)2 (2m+l)2 (2/1+1)2

This expression requires the summation of many terms at small values of time, but only a few terms at large values of time. The implementation of this complicated equation had to be checked against simpler

analytic cases. The first case represents a homogeneous solid occupying the semi-infinite region, x >

0. The solid is initially at unit temperature throughout. At time t = 0 the temperature at x = 0 is instantly reduced to zero. The temperature within the region is given by Carslaw and Jaeger (1959),

p. 59, equation (3): TM = erf (A9)

The second case involves a solid which occupies the region x > 0, y > 0, z > 0. It is initially at unit

temperature and at time t = 0 the temperature at the x = 0, y = 0, and z = 0 surfaces is instantly

reduced to zero. The temperature within the region is given by Carslaw and Jaeger (1959), p. 184,

equation (1):

T{x,y,z,t) = erf [ X \ erf [ y

] erf

r z \ (AlO)

The temperatures near the comers of the block should be very similar to this.

There was good agreement (error < 1 %) for a test with a = 30 mm, b = 20 mm, c = 10 nmi, and using

a 1 nmi uniform grid. The original variable grid model was found to be insufficiently accurate at points

where the grid size changed. It was therefore replaced by the current uniformly increasing grid at a cost of some increase in code complexity; although the results are not quite as accurate as for the uniformly

spaced grid, the diff'erence is very minor.

Test 5: One-Dimensional Transient Conduction, Two Different Materials

This test consists of one material in the region 1 (0 < z < L) initially at unit temperature and another

material in region 2 (z > L) initially at zero temperature. The boundary at z = 0 is adiabatic. At time

t = 0 heat begins to be conducted between the two regions. The temperatures in the two regions are given by Ozisik (1980), p. 328, equation (8-109):

A-4
where

11=0

(2«+l)I-:c

+ erfc

(2n+l)I+;c

n-O

2nI+ti(x-I)

- erfc

(2n+2)L+\xix-l)

(All)

There was good agreement between the FDA and analytic solutions. It even worked well when the first layer was only one-half of a grid thick. This may be useful for thin fabric coverings.

Test 6: Three-Dimensional Transient Conduction, Constant Heat Flux Impinging on a Circular

Area (Disk) on a Semi-Infinite Region.

The impinging flux is q, there is convective cooling from the surface, h(Jg — TJ, and the disk radius

is R. The temperature at the center of the heated area is given by Thomas (1957), equation (5):

< \\ TtR

4at

1 -exp

4at/J

+ erfc

- 2£ f

h J

—y/at

K

(A12)

«-o

1 - exp

[4kW)\

2 wtf" erfc((<>) Jo)

The first part of the above expression is the center point temperature if there is no convective cooling. Tests indicate that the accuracy of the FDA for this test is primarily dependent on how accurately the

circular flux pattern is represented on the rectangular surface grid. A small grid and assigning cell heat

gain according to the portion of the cell that is within the circle improve accuracy.

Note that tests 1 through 6 involve a step change, which should be the worst condition to simulate with

the FDA. In all cases the maximum errors occurred at the first time step, and the error declined as time

increased.

Test 7: Three-Dimensional Transient Conduction, Uniformly Moving Point-Source Heat Flux

This test consists of a point source of power Q moving in an infinite body at a constant velocity v in the

X direction (Schneider (1973), pp.3-86, equation (78)): —(T-TJ = -^exp

2a

(A13)

A-5
where ^ = x - vt and r = (^^ + + z^)''^. Adjusting for a semi-infinite body with the point source moving along an adiabatic surface is done by replacing 47r by Itt in equation (A 13). Generally good

agreement was achieved for this test. Accuracy was limited by grid size near the point source and the

fact the this is a quasi-steady case in that movement of the point source does not have a beginning point.

Larkin's Method

The FDA algorithm in CTEST3 was transferred directly into the new substrate model TMPSUB2. The

addition of pyrolysis forced the use of a very small grid in the region of peak temperature for a satisfactory solution. This, combined with the stability requirement of the explicit Euler method, forced a very

small time step and therefore a very long execution time. A different FDA algorithm had to be found

to achieve a program fast enough to be useful. Larkin's method was chosen because of its simplicity in

that it uses the same spatial FDA as the original method while the new temporal FDA does not require

the solution of simultaneous equations.

CTEST3 was not rewritten to run all of the test cases, but several tests were made with TMPSUB2 (some

using a modified surface boundary condition) which indicate the accuracy of the method for different grid

size and time step options. The results of these tests are shown in Tables A-1, A-2, and A-3.

Table A-1 gives the results of several tests which can be compared to equation (A5) to determine the

effect of grid spacing. These tests use a region 40 mm thick to simulate a semi-infinite body which is shown to be appropriate by having negligible heat flux at the constant temperature surface at Z=40 mm.

Comparisons are made based on the temperature of the surface (Z=0 mm).

Test 2a: Using a constant 1 mm grid spacing the temperature after 100 seconds is 0.33% less than the exact value.

Test 2b: Using a constant 0.5 mm spacing gives a surface temperature 0.18% below the

theoretical value.

Tests 2c through 2f use variable grid spacing to reduce the number of cells and execution time

at the cost of some loss of accuracy.

Table A-2 gives the results of several tests where the parameters that control the time step while

maintaining a constant grid spacing, are varied. These parameters are the maximum time step, dtmax,

and the maximum temperature change, dTmax. (Whenever T„+j - T„ exceeds dTmax, the time step is halved.) Obviously the greatest accuracy should be achieved with small values for these two parameters,

but execution time is reduced by using large values. There is no obvious optimum; the user must choose

values appropriate for results he wishes to achieve.

Table A-3 shows tests of different grid spacings for full three-dimensional heat conduction ft-om a

stationary spot heat flux. These tests use representative thermal properties for the fabric and padding.

Various combinations of cell spacings are used to select the best grid.

A-6
Table A-1. Grid Spacing Tests

Thermal diffusivity: 2e-07 v?

/

b Thermal conductivity: 0.1 W/mK

Surface heat flux: le+04 W/m^ Distance from surface: 0 m

time Texact Test2a Test2b

0000 0. 0000 0 . 000 0 .000

0125 178412 5 . 000 9 . 319

0250 252313 9 . 643 16 . 889

0. 375 30. 9019 13 . 970 23 . 331

0500 356825 18 . 298 28 .916

0

625 398942 22 .252 33 . 840

0750 43

7019 26 . 206 38 . 035

0

875 472035 29 . 838 42 .229

1000 504627 33 . 469 45 . 757

1250 564190 40 . 170 52 . 366

1500 618039 46 . 383 58 . 208

1. 750 66. 7558 52 . 167 63 . 496

2000 713650 57 . 577 68 . 361

2250 75 . 6940 62 . 656 73 . 042

2500 79 . 7885 67 . 445 77 . 173

2750 83 . 6828 71 . 859 81 . 305

3000 87 . 4039 76 . 273 85 . 049

3500 94

4070 84 . 279 92 . 243

4

000 1009253 91 . 625 98 .911

4500 107 . 0474 98 . 433 105 . 154

5 . 000 112 . 8379 104 .796 111 . 046

6

000 123

6077 116 . 456 121 . 976

7 . 000 1335115 127 . 018 132 . 004

8000 142 . 7299 136 . 550 141 .315

3 . 000 TCI 151

3880 145 .711 ICO 150 . 034

10 . 000 159 . C T £ Q

5 7 69 154 O O Q • z29 ICO 158 O T O• 272

15 . 000 19 5 . 4410 191 . 150 194 • 315

20. 000 22 5. ^ T tr o 57 58 221 n o ^

. 985 224 c ^ r\. 570

25. 000 4b «J ^ . ij M iJ £, ^ ^ 7 * \J -J ^ X • 7 U 30. 000 276. 3953 273 .438 275 .472

35. 000 298. 5410 295 . 758 297 .599

40. 000 319. 1538 316 .566 318 .155

45. 000 338. 5137 336 .042 337 .525

50. 000 356. 8248 354 .491 355 .815

60. 000 390. 8820 388 .738 389 .887

70. 000 422. 2008 420 .206 421 .230

80. 000 451. 3517 449 .479 450 .409

90. 000 478. 7307 476 .960 477 .816

100. 000 504. 6265 502 .943 503 .740

cells

41 81

steps 153 154

Test2c Test2d Test2e Test2f

0 .000 0 .000 0 .000 0 .000

9 .319 9 . 319 9 . 319 9 .319

16 .889 16 .889 16 .889 16 .889

23 .331 23 .331 23 .331 23 .331

28 .916 28 .916 28 .916 28 .916

33 .840 33 .840 33 .840 33 .840

38 .035 38 .035 38 .034 38 .034

42 .229 42 .229 42 .229 42 .229

45 .757 45 .757 45 .757 45 .757

52 .366 52 . 366 52 .366 52 .366

58 .208 58 .208 58 .208 58 .208

63 .496 63 .496 63 .496 63 .496

68 . 361 68 .361 68 .361 68 .361

73 .042 73 .041 73 .041 73 .041

77 . 173 77 . 173 77 . 173 77 . 171

81 .305 81 .304 81 .304 81 .302

85 . 049 85 .048 85 . 047 85 .044

92 . 243 92 . 242 92 . 240 92 .233

98 .909 98 . 908 98 .904 98 . 893

105 . 152 105 . 149 105 . 143 105 . 125

111 . 042 111 . 038 111 .030 111 .002

121 . 970 121 .962 121 . 946 121 . 894

131 . 993 131 . 981 132 . 050 131 .966

141 . 300 141 . 283 141 . 309 141 . 188

150 . 014 149 . 991 149 . 993 149 . 830

158 .248 158 .219 158 .201 157 . 994

194 . 266 194 . 208 194 . 109 193 . 654

224 . 600 224 . 514 224 . 343 223 . 631

251 . 307 251 . 195 250 .960 249 . 998

275 .367 275 .235 274 .852 273 .660

297 .487 297 .339 296 .828 295 .425

318 .040 317 .879 317 .378 315 .771

337 .408 337 .234 336 .650 334 .844

355 .694 355 .509 354 .928 352 .929

389 .758 389 .550 388 .899 386 .532

421 .093 420 .862 420 .147 417 .431

450 .262 450 .010 449 .235 446 .188

477 .660 477 .388 476 .556 473 .195

503 .574 503 .282 502 .397 498 .738

35 25 18 12

154 154 153 153

TEST2A: dz=1.0min,

TEST2B: dz=0.5mm,

TEST2C: dz=0.5mm,

TEST2D: dz=0.5mm,

TEST2E: dz=0.5mm,

TEST2F: dz=0.5mm,

nz=41, nc=41,

nz=81, nc=81,

nz=35, nc=4,

nz=25, nc=4,

nz=18, nc=4,

nz=12, nc=4,

r=1.000, dtmax=l

r=1.000, dtmax=l

r=1.055, dtmax=l

r=1.115, dtmax=l

r=1.234, dtmax=l

r=1.628, dtmax=l

.0, dTmax = 5.0

.0, dTmax = 5.0

.0, dTmax = 5.0

.0, dTmax = 5.0

.0, dTmax = 5.0

.0, dTmax = 5.0

A-7
Table A-2. Time Step Control Tests

Thermal diffusivity: 2e-07 m^/s Thermal conductivity: 0.1 W/mK

Surface heat flux: le+04 W/m^ Distance from surface: 0 m

time Texact Test2b Test2g Test2h Test2J

0. 000 0. 0000 0 . 000 0 . 000 0 .000 0. 000

0. 125 17. 8412 9 .319 9 . 093

0. 250 25. 2313 16 .889 16 . 685 20 . 000 20. 000

0. 375 30. 9019 23 .331 23 . 148

0. 500 35. 6825 28 .916 28 .750 31 .538 31. 538

0. 625 39. 8942 33 .840 33 .688

0. 750 43. 7019 38 .035 38 . 107 39 . 690 39. 690

0. 875 47. 2035 42 .229 42 . 125

1. 000 50. 4627 45 .757 45 .801 47 .841 47. 841

1. 250 56. 4190 52 . 366 52 .390

1. 500 61. 8039 58 .208 58 .223 58 .553 58. 553

1. 750 66. 7558 63 .496 63 .506

2 . 000 71. 3650 68 .361 68 .368 69 .264 69

.

264

2. 250 75. 6940 73 .042 72 .899

2 . 500 79. 7885 77 . 173 77 . 158 77 .403 77. 403

2 . 750 83 . 6828 81 . 305 81 . 191

3. 000 87 . 4039 85 .049 85 . 028 85 . 542 85. 542

3. 500 94. 4070 92 .243 92 .222

4. 000 100. 9253 98 .911 98 .891 97 .831 97. 831

4 . 500 107. 0474 105 . 154 105 . 136

5. 000 112. 8379 111 .046 111 .029 110 . 120 110. 120

6. 000 123. 6077 121 .976 121 .963 120 .390 120. 390

7. 000 133 . 5116 132 .004 131 .993 130 .660 130. 660

8. 000 142. 7299 141 .315 141 .312 139 . 622 139. 622

9. 000 151. 3880 150 .034 150 .053 148 .583

10. 000 159. 5769 158 .272 158 .312 156 . 621 157. 545

15. 000 195. 4410 194 .315 194 .402 192 .945

20. 000 225. 6758 224 .670 224 .773 223 .337 221. 615

25. 000 252. 3132 251 .396 251 .504 250 .229

30. 000 276. 3953 275 .472 275 .656 274 .420

35. 000 298. 5410 297 .599 297 .855 296 .726

40. 000 319. 1538 318 .155 318 .512 317 .415 313. 420

45. 000 338. 5137 337 .525 337 .900 336 .887

50. 000 356. 8248 355 .815 356 .236 355 .255

60. 000 390. 8820 389 .887 390 .336 389 .440 382 . 464

70. 000 422. 2008 421 .230 421 .690 420 .860

o •

00 000 451. 3517 450 .409 450 .870 450 .093 442 . 530

90. 000 478. 7307 477 .816 478 .274 477 .541

100. 000 504. 6265 503 .740 504 .191 503 .495 495. 723

cells: 81 81 81 81

steps: 154 724 106 42

TEST2B: dz=0.5mm, nz=81, nc=81, r=1.000, dtmax=1.0, dTmax = 5.0

TEST2G: dz=0.5mm, nz=81, nc=81, r=1.000, dtmax=1.0, dTmax = 1.0

TEST2H: dz=0.5mm, nz=81, nc=81, r=1.000, dtmax=1.0, dTmax = 20.0

TEST2I: dz=0.5mm, nz=81, nc=81, r=1.000, dtmax=4.0, dTmax =5.0

( Same results as test2b because dT > 2.5 at dt = 1.0 ) TEST2J: dz=0.5mm, nz=81, nc=81, r=1.000, dtmax=4.0, dTmax = 20.0

A-8
TaDle A- 3. 3-D Transient Conduction Tests

time test3a test3b test3c test3d test3e

0. 00 0.000 0 .000 0 . 000 0.000 0.000

1. 00 65.749 70 .904 71 . 046 71.254 71. 301

2. 00 88.450 91 . 190 91.381 91.660 91.723

3. 00 104 . 807 106 .828 107 . 053 107.382 107.455

4. 00 118.53 6 120 .437 120 . 692 121.060 121. 140

5. 00 13 0.745 132 . 693 132 . 974 133.375 133 .460

10. 00 178 .2 69 180 .414 180 .785 181.287 181. 377

15. 00 1 T AAA 211 .444 213 .792 214.214 214.765 214.850

20. 00 2 3 6. 140 238 . 668 239 . 121 239.692 239.766

3 0 . 00 2 59 . 7 52 272 . 068 272 . 539 273 . 119 273 . 173

40. 00 291 . 304 293 .512 293 . 988 294.550 294.585

50. 00 305 .937 308 . 182 o /^ o f c c\ 308 . 659 309 . 201 309 . 219

60. 00 316.279 318 . 588 319. 064 319.589 319.594

70. 00 323 .830 326 . 195 326.670 327.180 327. 176

80. 00 329.498 331 . 903 332.376 332.876 332.865

90. 00 o o o o c o 333od3 336 .285 335 . 757 337.248 337,233

100. 00 337.271 339 .718 340 . 188 340.673 340. 654

cells18081 18081 18081 19176 18375

steps152 154 154 154 154

time • 303 . 08 316 .37 316. 10 338.90 346. 32

TEST3A: dx=0.5mm, nx=41, ny=21, nz = 21, nc = 4, yw=zw = 40,

rx=ry=1.167, rz=1.187, dtmax=1.0, dTmax = 5.0

TEST3B: dx=0.25mm, nx=41, ny=21, nz = 21, nc = 4, yw=zw = 40,

rx=ry=1.240, rz=1.240, dtmax=1.0, dTmax = 5.0

TEST3C: dx=0.25mm, nx=41, ny=21, nz = 21, nc = 4, yw=zw = 30,

rx=ry=1.210, rz=1.210, dtmax=1.0, dTmax = 5.0

TEST3D: dx=0.25mm, nx=47, ny=24, nz = 17, nc = 4, yw=zw = 30,

rx=ry=1.161, rz=1.32 6, dtmax=1.0, dTmax = 5.0

TEST3E: dx=0.25mm, nx=49, ny=25, nz = 15, nc = 4, yw=zw = 30,

rx=ry=1.149, rz=1.431, dtmax=1.0, dTmax = 5.0

A-9

APPENDIX B

TMPSUB2 AND SUBSTRAT USER'S GUIDE

General Information

SUBSTRAT and two helper programs are

available on a diskette for IBM PC compatible computers. The distribution diskette includes both executable and source code for SUBSTRAT, executable code for

the other two programs, and sample input

files. The general relationship of programs and files is illustrated in Figure B1. The TSDATA program is used to prepare data files for SUBSTRAT which

in turn creates two types of output files. The plot file is used by the TSPLOT

program to display contour plots of substrate temperatures. The list file includes

a step-by-step record of the highest temperature in the substrate. Figure B-1. TMPSUB2 Programs and

TSDATA and TSPLOT must be run using Files

MS-DOS on IBM PC compatible computers with VGA graphics. SUBSTRAT requires a 386 class PC with math coprocessor or 486 class PC

(no graphics needed) in order to achieve satisfactory performance. Typical execution times are 20 to 30

minutes on a 33Mhz 486 computer.

The SUBSTRAT source code (file SUBSTRAT.CCC) can be compiled using any ANSI C compiler. All

SUBSTRAT input and output files are ASCII files. Therefore, SUBSTRAT can be recompiled and run on a diff'erent computer, while still using TSDATA and TSPLOT on a PC and transferring files between

the computers. A diff"erent computer may allow SUBSTRAT to execute faster and/or handle more cells for more accurate simulation.

Be sure to inspect the README file on the distribution diskette. One way to read this file is to place

the diskette in drive A: (or drive B:) and type

MORE < A:README (or MORE < B:README)

A permanent copy may be made with

PRINT <A:README

The README file contains a list of all files on the diskette, instructions for installing the necessary files on your hard disk, and information on any changes or additions to the program.

There should be at least 1,000,000 bytes available on your hard disk. It is best to create a single

subdirectory for the executable programs and related data files. This will allow you to easily delete all

TSDATA

I

data

TMPSUB2

list plot

::TS:Biii::

B-1
the files related to this program when you are finished with it. In general, when running on a PC, keep

all files in the current working directory.

The following sections give details of the operation of SUBSTRAT, TSDATA, and TSPLOT.

SUBSTRAT

All input for SUBSTRAT must be placed in a data file. This file is read as the "standard input stream,"

so SUBSTRAT executes on MS-DOS and UNIX computers by redirecting the input file. For example,

if there already exists an input data file, e.g., TEST1.DAT, one can begin a run by typing

SUBSTRAT < TEST1.DAT

Note: SUBSTRAT can be aborted at any time by simultaneously pressing CTRL and C.

Sample run

If you have not already done so, install the SUBSTRAT program following the instructions in the

README file. For example, from the directory on the hard disk where you want SUBSTRAT installed

and with the diskette in drive A:, type

A:INSTALL A:

At the present time you only need to install in response to the first question you will be asked.

During installation a program will determine the "initial unallocated memory." If the required memory

is greater than this amount, follow the instructions under memory requirements below before proceeding.

Should you wish to examine the directory at this point, a DIR command should indicate the presence of

at least the following files: SUBSTRAT.EXE, TEST1.DAT, TEST2.DAT, and MEMREM.EXE. As

mentioned above, begin the sample run by typing

SUBSTRAT < TEST1.DAT

SUBSTRAT will then (1) indicate that it is reading the data file, (2) echo the simulation title to the

screen, (3) initialize the data in all cells, and (4) report the amount of unallocated memory. If there is "insufficient memory," an error message will be displayed and the run aborted. See the memory

requirements section below for corrective action.

Every time SUBSTRAT completes a time step during the simulation, it displays the time (in seconds) and

peak temperature (degrees C) on the screen. This allows you to monitor the progress of the simulation.

This same information is included on the list file, TESTl.LST, thus providing a permanent record of the

primary value computed during the simulation. The list file also includes an echo of the input file which

helps to identify the simulation and is especially useful in identifying any errors in the input which are

also written to the list file.

B-2
The ploi file, TESTl.PLT, contains all temperatures in the Y=0 and Z=0 planes at the times specified

plus data to identify the simulation and the coordinate values.

Memory Requirements

SUBSTRAT has been written to handle an arbitrary number of cells up to some limit imposed by

available memory or by the operating system. MS-DOS limits available random access memory (RAM,

not disk memory) to 640,000 bytes which is shared by the program, the data in the program, a portion

of the operating system, and perhaps various TSR (Terminate & Stay Resident) programs. TMPSUB

allocates memory for its data arrays from available RAM. When there is insufficient memory to run a

particular simulation, you must either make more memory available or reduce the memory required by

the simulation. One way to make more memory available is to remove TSR programs, such as network

connections. This usually requires rebooting the computer in such a maimer that these programs are not

automatically accessed. MS-DOS 5.0 uses somewhat less memory than previous versions. See your PC

consultant for assistance. Run MEMREM to determine approximately the memory available for data.

The memory required by the simulation is determined primarily by the total number of cells which in turn

is determined by the values on line 4 of the data file (see below). The total number of cells can be

reduced to an arbitrarily small value, but accuracy will sufl'er. If you must search for the number of cells than can be run on your machine, begin with a relatively small number and increase it until the

unallocated memory reported by SUBSTRAT is small.

With minimal other uses of memory, SUBSTRAT is limited to less than 20,000 cells under MS-DOS.

The only way to simulate more cells is to recompile SUBSTRAT for a diff'erent operating system and/or

a different computer.

B-3
Contents of SUBSTRAT Data File

line variables brief description

1 LIST Name of the output file

2 PLOT Name of the plot file

3 TITLE Project title; echoed to output

4 NX NY NZ number of cells in the X, Y, and Z directions

( total number of cells = NX x NY x NZ )

5 NCX DXO XW XO X-coordinate data':

NCX number of constant length nodes

in both directions from XO

DXO length [mm] of constant length nodes

XW total length [mm] of X axis

XO location of center of region of constant length nodes

= initial position of peak flux from cigarette

6 NCY DYO YW Y-coordinate data':

NCY number of constant width nodes from Y = 0

DYO width [mm] of constant width nodes

YW total width [mm] of Y axis

7 NCZ DZO ZW FT Z-coordinate data':

NCZ number of constant depth nodes from Z = 0

DZO width [mm] of constant depth nodes

ZW total depth [mm] of Z axis

FT fabric thickness [mm]^

8 SEP Width of air gap between fabric and padding [mm]

9 E TM NCF Fabric data:

E emittance

TM maximum temperature^ [°C]

NCF number of data points for T, K, C (NCF <= 10)

10 T K C Fabric thermal properties*:

T temperature ["C]

K conductivity [W/mK]

C specific heat [kJ/kgK]

Repeat line 10 up to NCF times.

11 DV DC DA Densities:

DV virgin material [kg/m^]

DC char [kg/m^]

DA ash (kg/m-']

12 A Nl N2 T H Fabric non-oxidative pyrolysis data:

A reaction rate coefficient [1/min]

Nl fabric mass exponent

N2 oxygen concentration exponent (must be zero)

T activation temperature [K]

H heat of pyrolysis [kJ/kg]

13 A Nl N2 T H Fabric oxidative pyrolysis data:

A reaction rate coefficient [1/min]

Nl fabric mass exponent

N2 oxygen concentration exponent

T activation temperature [K]

H heat of pyrolysis [kJ/kg]

B-4
14 A Nl N2 T H Fabric char pyrolysis data (similar to 13)

15 E TM NCP Padding data:

E emittance

TM maximum temperature^ ["C]

NCP number of data points for T, K, C (NCP <= 10)

16 T K C Padding thermal properties^:

T temperature [°C]

K conductivity [W/mK]

C specific heat [kJ/kgK]

Repeat line 16 up to NCP times.

17 DV DC DA Densities:

DV virgin material [kg/m^]

DC char [kg/m^]

DA ash [kg/m^]

18 A Nl N2 T H Padding non-oxidative pyrolysis data (similar to 12)

19 A Nl N2 T H Padding oxidative pyrolysis data (similar to 13)

20 A Nl N2 T H Padding char pyrolysis data (similar to 13)

21 BC TA TO OX Boundary & initial conditions:

BC 0 = adiabatic outer boundaries^,

1 = constant temperature outer boundaries

TA ambient temperature [°C]

TO initial substrate temperature [°C]

OX oxygen concentration [fraction]

22 QO V Y X+ X- Moving radiant flux distribution on surface:

QO peak flux [kW/m^]

V velocity [mm/min]

Y standard deviation in Y direction [mm]

X+ std dev in positive X direction [mm]

X- std dev in negative X direction [mm]

23 HC HG Heat transfer coefficients^

HC for quiescent air [W/m^K]

HG for impinging air [W/m^K]

24 dt DT TT Simulation control:

dt maximum time step [s]

DT maximum temperature change [°C]

TT total simulation time [s]

25 LIST I J K Debug reports (SUBSTRAT compiled with DEBUG defined):

LIST 1 = activate data dumps (default 0)

I J K are the indices of the cell to be studied

26 PLOT Time at which to write temperatures to the plot file;

repeat line 26 up to 10 times.

Notes:

' The relationship between the number of nodes, number of constant size nodes, size of constant size

nodes and total length of an axis must be such that the variable width nodes are increasing in size.

B-5
^ The fabric thickness and the depth of the constant depth cells are related: the fabric/padding boundary

must be halfway between two cells. Therefore, the fabric thickness must be 0.5, or 1.5, or 2.5, etc., times the depth of the constant depth cells, or the constant depth cells must be 2, or 2/3, or 2/5, etc. times the thickness of the fabric. The boundary between the fabric and padding must be within the region

of constant depth cells.

^ When the temperature of any fabric (or padding) cell reaches the "maximum" fabric (or padding)

temperature, the simulation is terminated.

The way that k(T) and c(T) are input is by entering the values for each at a number of temperatures;

the program then carries out a cubic spline fit to those points, to obtain the values at any other

temperature. Even if k(T) and c(T) are given by explicit equations, this is still the way that the program

"knows" the values.

^ The "outer boundaries" of the substrate consist of the X = 0, X = XW, Y = YW, and Z = ZW

planes (see items 5, 6, & 7).

^ These coefficients refer to cases A and B in Section II. B. 6. When using HC set HG to zero and vice

versa.

The contents of the data file are described further in the description of the TSDATA program. It will

generally be easier to check a data file by processing it with TSDATA, than to compare it line by line

with the description given above.

B-6
TSDATA

The SUBSTRAT data file can be created with any ASCII line editor. However, the contents of the file are quite cryptic and thus prone to error; far better is to use TSDATA. TSDATA is an interactive

program for creating data files. Because it is interactive, it uses certain commands which restrict its operation to IBM PC compatible computers. It includes extensive checking of the input data. TSDATA

is especially useful for creating a data file which is only slightly diff'erent from another data file. This

is useful in performing the parametric studies for which SUBSTRAT was designed.

Two special files are used by TSDATA. The help file, TSDATA.HLP, contains the text of the

interactive help messages. Help is activated by pressing the Fl function key. If the help file is not

available in the current working directory, no interactive help will be available. The configuration file, TSDATA.CFG, sets the colors of the display. The file included on the distribution diskette assumes that a standard VGA monitor is being used. If the configuration file is not in the current working directory, a set of default colors will be used. A new configuration file can be made by using the MAKECFGT

program. See the README file for instructions.

The operation of TSDATA is explained on the following pages which show the messages and input screens which will appear as the program is run. After reading through these pages, try using TSDATA

with one of the sample data files. Begin the program by typing TSDATA. Abort the program by

pressing CTRL and C.

Sample Runs

At this point the following files should be available in the current working directory: TSDATA.EXE,

TSDATA.HLP, TSDATA.CFG, and TEST1.DAT. If they are not, use the INSTALL procedure.

First use TSDATA to view the contents of the TEST1.DAT file. Type:

TSDATA

Press ENTER until the main menu appears; press ENTER again to set up the file information. Enter

TEST1.DAT as the name of the previous data file. Press ENTER to cycle through file names (respond

Y to the warning message about a duplicate file). Press ESC to return to the main menu. Press ENTER

to view each section of data in turn; press ESC to return. When you reach the Save option, press CTRL

and C to abort the program. This process will not have changed TEST1.DAT.

Now create a new data file, TEST3.DAT, similar to TEST1.DAT but with a peak flux of lOkW/m^.

Proceed as above until you reach "Name of new data file:". At this point move the cursor to the T in

the file name and press 3 and then press ENTER. Revise the other two file names and the title accordingly. Press ESC to return to the main menu, and move to the boundary conditions option. Press ENTER and move (using cursor keys) to the peak heat flux. Change the value to 20, press ENTER,

press ESC. Now enter the save option, and respond Y to save the file TEST3.DAT. Then enter the exit

option, respond Y to exit, and respond Y to create a batch file. After this you should have exited

TSDATA. Now type RUN to start SUBSTRAT using TEST3.DAT as input data.

B-7
TSDATA Input Screens and Help Messages

Standard title page with disclaimer at start of TSDATA program.

TSDATA — interactive program to prepare data files for SUBSTRAT

Version 1.0

Developed at the National Institute of Standards and Technology.

Progreun author: George Walton

This program is furnished by the government and is accepted by any

recipient with the express understanding that the United States

Government makes no warranty, expressed or implied, concerning the

accuracy, completeness, reliability, usability, or suitability for

any particular purpose of the information and data contained in

this program or furnished in connection therewith, and the United

States shall be under no liability whatsoever to any person by reason

of any use made thereof. This program belongs to the government.

Therefore, the recipient further agrees not to assert any proprietary

rights therein or to represent this program to anyone as other than

a government program.

General description of input process:

This program assists you in preparing input data files for the SUBSTRAT

progreim. It operates best by reading an existing SUBSTRAT data file

which is then modified to create a new data file. A sample data file

is distributed with the program. It can also be used to enter data

from scratch. Several data files can be created in one TSDATA session.

Data are processed interactively through a system of data entry menus. Keyboard input is required from the user whenever the cursor is in

a data entry field. A data entry field is designated by a special

color as is shown in the lower right corner of this screen. This is

a standard pause allowing the user to read the screen. Pressing any

key in response will allow the program to continue.

While the cursor is in a data entry field, it will often be possible

to get help by pressing the Fl function key. Help is intended to

give additional information about the data. For help to work, the

TSDATA. HLP file must be in the same directory as the TSDATA.EXE file.

Progreun execution may be terminated when the cursor is in a data

entry field by pressing CTRL and C simultaneously.

B-8
General description (continued):

There will usually be several data entry fields on a single screen. The field may be blank or it may present a default response. The

contents of a data entry field (even all blank) may be edited.

That is, characters may be overwritten, deleted, or inserted.

It may or may not be necessary to make an entry depending on context.

These fields are ordered from top to bottom and left to right.

Pressing the tab key or the down-arrow key moves the cursor to the

next field to the right/down; The shifted-tab or up-arrow moves

to the previous field to left /up. Control-home or page-up moves the cursor to the first field on the screen. Control-end or page-down

moves to the last field. Movement between fields can be done only

if an entry is not required and no other keys have been pressed.

Data entry begins in the exchange mode, i.e., the value of the key

pressed replaces the character at the cursor. Pressing the insert

key will switch to the insert mode (and from insert to exchange).

The delete key will remove the character at the cursor. Control-x

will clear the entire data entry field. Move the cursor left with

the left-arrow, control-left-arrow, or home keys. Move the cursor

right with the right-arrow, control-right-arrow, or end keys. When

the data is satisfactory, press the ENTER key.

Various checks are usually made on each data entry. These may produce

a warning or error message at the bottom of the screen. The first line

of this message indicates the nature of the problem. The second line

indicates the severity of the problem, the file and line in the TSDATA

source code where the error message originated, and whether some help

may be available by pressing the Fl key.

Most errors will return you to the data entry field to correct the input.

Such errors include invalid characters or numeric values outside certain

situation-dependent limits.

Some problems may generate a question:

Question? (y/n) ^

The user can press the Y and then ENTER keys to indicate a positive (yes)

response or N and then ENTER to indicate a negative (no) response.

Some errors are fatal causing the program to terminate.

B-9
Screen 1: (primary menu)

SUBSTRAT data preparation: ^ File information ^ Geometry description ^ Fabric data ^ Padding data ^ Boundary conditions ^ Simulation control ^ Save this data file ^ Exit data preparation

Use cursor keys to move between menu selections.

Press ENTER to activate the menu selection at the X.

Press ESC to return from a selection. Press Fl for help.

General help message:

This program assists you in preparing input data files for the SUBSTRAT

program.

This is the main menu. It directs you to the different data

preparation subsections.

Begin each data file by entering the "file information".

Data relating to the geometry, fabric, padding, or boundary conditions

may be entered or changed in any order.

The data file is not created until you "save this data file".

Multiple data files can be created in one interactive session.

Terminate the session by entering "exit data preparation".

Note:

An X appears in the selection designator: |.

B-10
Screen 2:

File information: press ESC when done; press Fl for help.

Copy data from previous data file ncimed:

Name of new data file: IMiliii^^

Name of list file:

Name of plot file:

Title for this simulation:

Files previously defined in this session:

data files list files plot files titles (first 40 characters)

I3:i8iitt»iuiaiiiiu:!:ii:::!>ii !::i:!;::ii:i:!:::t:!ii:i:i!i!:;:iii ;:i»!:::!:!:n:::i;:mmm>m» n!:ii::iiiii:iiiiti:iiiiiii:iii:iiiti!:!:it!:!!::::::

!S!n!2!H!H!!l!!!!!!! itiii::!:ii'.iiii !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

iiiiiii!:iiiii:tit!ittiist:iiii!:

•.!:!::::!:!:!ti!:!:!! i::::!:::!!:!!::ti

!!i:in!!i:iiii!:!

General help message:

The data files you create are stored in the current directory. If you

want to make a data file which is similar to a previous data file, enter

the name of the previous file. That file must also be in the current

local directory.

SUBSTRAT creates two files when it is run:

(1) the list file (which echoes the input, saves any error meesages, and

notes the peak temperature of the substrate during the simulation)

,

(2) the plot file (which becomes the input file to the plotting program

used to view the temperature distribution in the substrate).

Files should generally have different names. For exsumple, consecutive

SUBSTRAT runs with the same plot file names will save only the plot file

from the last run (which will have replaced the previous plot files).

Use identical names if you definitely want to replace existing files,

including data files.

The title that you enter here is echoed in the list and plot files. Use

the title as a reminder of the special features of a particular run.

Up to 10 data files created previously in this session are recorded at

the bottom of the screen for your convenience.

B-11
Help messages in response to error messages:

(1) could not open file:

The file which you specified was not found. A typing error is the most

likely problem, or the file is not in the current directory.

(2) duplicate file names:

This message indicates that a file name matches one in the local directory,

or one previously set in this session, or that you have not pressed ENTER

at each of the file names to complete the check of file names.

If you have copied an existing data file, the "new" file neunes are the ones that appeared in the existing data file.

Files should generally have different naunes. For example, consecutive

SUBSTRAT runs with the Seime plot file names will save only the plot file

from the last run (which will have replaced the previous plot files).

Use identical names if you definitely want to replace existing files,

including data files.

You may use a duplicate file neune by responding 'Y' to the question

Gibout writing over the previous file.

Notes:

An entry for the previous data file causes its data to be copied.

The new data file will be created when "save this data file" is executed at the main menu. This save adds

the new data, list and plot file names to the list displayed at the bottom of this screen.

There should be no duplicate file names. That would cause files to be overwritten during simulation.

The program checks for duplicates.

The simulation title is echoed in the output files.

Note that upon entering the "Title for this simulation," the cursor jumps back to the first box and blanks

it out. This is done only so that if the user has changed his mind (or an error has been made), new (or

corrected) entries can be made immediately. If there has been no error, then simply press ESC.

B-12
Screen 3:

Geometric description:

Fabric thickness:

Width of air gap:

press ESC when done; press Fl for help. mmmm

X-coordinate data:

Number of grid points;

Total length of substrate:

Size of constant length cells;

Number of constant length cells:

Y-coordinate data:

Number of grid points!

Half-width of substrate!

Size of constant width cells!

Number of constant width cells!

Z-coordinate data:

Number of grid points!

Total depth of substrate!

Size of constant depth cells!

Number of constant depth cells:

Rx = 1 . XXX

Ry = l.xxx mm

mm

mm

mm

Rz = 1 . XXX

cells: xxxxx

General help message:

The substrate consists of a thin fabric, padding, and possibly an air gap between them. An air gap width of 0.0 indicates no gap.

The X direction is along the cigarette; the Y direction is along the

fabric; the Z direction is into the padding.

A variable grid is used. It consists of several constant width cells near the point of peak incident heat flux followed by increasingly

larger cells out to the boundaries of the substrate. For example:

I I I I I I 1 I I I I

I I I I I I I I I I I Y=0 Y=width

Each variable width cell is R times longer than the preceding cell.

Simulation is most accurate when R is only slightly greater than one. A value of R less than about 1.25 should be sufficiently accurate.

Simulation will generally involve a trade-off of accuracy and run time.

The fabric thickness and the depth of the constant depth cells are related: the fabric/padding boundary must be halfway between two cells.

Therefore, the fabric thickness must be 0.5, or 1.5, or 2.5, etc., times

the depth of the constant depth cells.

Other help messages:

For the X-axis the following conditions should apply:

4 * NCX < NX and (NX - 1) * DXO <= XW

where NX = number of grid points, XW = total length of substrate,

DXO = size of const length cells, & NCX = number of const length cells.

For the Y-axis the following conditions should apply:

B-13
2 * NCY < NY and (NY - 1) * DYO <= YW

where NY = number of grid points, 2YW = total width of substrate,

DYO = size of const width cells, & NCY = number of const width cells.

Because of the assumed bilateral symmetry of the flux, it is only necesary

to make calculations between y = 0 and y = YW (half the width). The number

of constant-width cells is that in the half-width section.

For the Z-axis the following conditions should apply:

2 * NCZ < NZ and (NZ - 1) * DZO <= ZW

where NZ = number of grid points, ZW = total depth of substrate,

DZO = size of const depth cells, & NCZ = number of const depth cells.

The fabric thickness and the depth of the constant depth cells are

related: the fabric/padding boundary must be halfway between two cells.

Therefore, the fabric thickness must be 0.5 or 1.5 or 2.5, etc., times

the depth of the constant depth cells, or the constant depth cells must

be 2 or 2/3 or 2/5, etc. times the thickness of the fabric.

The boundary between the fabric and padding must be within the region

of constant depth cells.

Notes:

"Rx", "Ry", and "Rz" are the geometric progression rates. The sample problems have a relatively high

value in tiie Z direction. Tests have indicated this is satisfactory because the low conductivity of the

padding permits less heat transfer in this direction.

"cells:" gives the total number of cells used to model the substrate. It equals NX x NY x NZ. This

value is critical to the total memory required for a simulation.

The following figure illustrates the substrate coordinate system and the variable grid:

Figure B-2 . Substrate Coordinate System and Variable Grid

B-14
Screen 4:

Fabric Properties: press ESC when done; press Fl for help.

Reactions

* A nl n2 Ta He

General help message:

Emittance is used in computing radiant heat loss from the fabric to

ambient and radiant heat transfer across the air gap.

The simulation stops when the temperature of any fabric cell reaches

the prescribed maximum temperature.

Thermal conductivity (K) and heat capacity (C) can be specified at up

to 10 temperatures. SUBSTRAT uses a cubic spline curve fit for K and C.

As fabric mass is lost by pyrolysis, K is also reduced in proportion

to the current density. Values are required at least two temperatures

even if K and C are constant.

A two-stage, three-reaction pyrolysis model is used: the virgin

material is converted to char in stage 1 and then to ash. The char

density is the value for completely converting the virgin material

to char with no conversion to ash. The ash density is what's left

after total pyrolysis.

The pyrolysis equations are of the form:

R = A * Dm'-nl * Cx"n2 * exp( Ta/T ) and Q = R * He

where

R = rate of pyrolysis

A = reaction coefficient [l/min]

Dm = density of material (virgin or char ) nl = related exponent

Cx = oxygen concentration

n2 = related exponent

Ta = activation temperature [K]

T = current cell temperature [K]

He = heat of pyrolysis [kJ/kg]

B-15
The first reaction is thermal degradation (n2 = 0 and He < 0) of virg

material. The second reaction is oxidation of the virgin material.

These two reactions produce char. The third reaction is oxidation

of the char to produce ash.

Help message in response to error message:

Temperatures must be given in increasing order.

The lowest temperature should be several degrees less than any

possible substrate temperature. The highest temperature should

be several degrees higher than the maximum temperature.

B-16
Screen 5:

Padding help messages are similar to fabric messages.

B-17
Screen 6: Boundary Conditions: press ESC when done; press Fl for help.

At X='0, X=Xmax, Y=Ymax, Z=Zmax: (M) adiabatic

(W) constant temperature

Initial substrate temperature:

Ambient temperature:

Oxygen mass fraction:

Heat transfer coefficients:

Quiescent air:

Impinging air:

Moving heat flux pattern from the^^^^i^arettej

Peak heat flux:

Initial X position of peak:

+X velocity:

±Y standard deviation (A)

:

+X standard deviation (B)

-X standard deviation (C) :

'C

'C

W/m^K

W/m^K

General help message:

Select the condition of the outer boundaries (X=0, X=length, Y=width,

Z=depth) by pressing ENTER in the appropriate place.

The incident heat flux from the cigarette is represented by a moving

flux distribution on the surface of the substrate. This distribution

has the shape of non-symmetric Gaussian curve which is described by six

parameters

:

(1) the X coordinate at the peak of the curve at the start of the

simulation { the initial position of the peak is (Xo, 0, 0) },

(2) the speed, S, at which the peak moves along the X axis { the

position P of the peak at time t is (Xo+S*t, 0, 0) },

(3) the maximum heat flux at the peak { at position (P, 0, 0) },

(4) the width of the curve in the Y direction { at position (P, ±A, 0),

the heat flux is 0.37 times the flux at the peak },

(5) the width of the curve in front of the peak { at position (P+B, 0, 0),

the heat flux is 0.37 times the flux at the peak }, and

(6) the width of the curve behind the peak { at position (P-C, 0, 0),

the heat flux is 0.37 times the flux at the peak }.

B-18
Screen 7:

Simulation Control: press ESC when done; press Fl for help.

General help message:

SUBSTRAT uses a variable time step which is chosen so as not to exceed

the maximum time step given and so that the maximum temperature change

in any cell is not greater than the maximum given. This allows the

program to run quickly when the temperatures are not changing rapidly.

Smaller values for these two parameters will lead to more accurate

simulations at the cost of longer execution times.

The simulation will stop at the total simulation time unless it has

already stopped by exceeding a maximum fabric or padding temperature.

Plot times are entered in increasing order. When the simulation reaches

a plot time, substrate temperatures are copied to the plot file. Leave

later positions blank if you do not want to use all 10 times. A plot is

automatically written when the simulation stops. The last plot time

should be greater than the simulation time to satisfy a requirement in

the SUBSTRAT program.

B-19
Screen 8:

Save this data file? (y/n) ^

display any error message (s).

General help message:

If you entered this area accidentally and are not ready to stop

preparing the current data file, respond "N" to this question.

This will return you to the main menu.

Help messages in response to error message:

(1) check data.

Data in the indicated section may be incomplete or incorrect.

Enter that section and check the data.

(2) temperature problem.

The lowest temperature for fabric or padding thermal properties must

be several degrees lower than either the initial or ambient temperature

You need to change (at least) one of those values.

(3) files limit.

No more SUBSTRAT data files can be created in this session.

You probably should exit TSDATA now. You may replace one of the

data files already created.

"i

B-20
Screen 9:

Exit data preparation? (y/n) ^ (1)

display error messages, if applicable

prepare DOS batch file to run all cases (y/n)? ||

(2)

(3)

Help messages:

(1) initial message.

If you entered this area accidentally and are not ready to stop

preparing data files, answer "N" to this question.

(2) error messages.

This warning indicates that TSDATA may have data which you have not saved.

If you respond "Y" to the question about checking this data file, you

are returned to the main menu where you can check and then save the file.

If you respond "N", you will continue to exit the TSDATA program.

(3) batch file message.

A positive response will create a file RUN. BAT resembling:

SUBSTRAT <datafile.l

SUBSTRAT <datafile.2

SUBSTRAT <datafiie.3

SUBSTRAT <datafile.4

This file can be used to run SUBSTRAT on computers using MS-DOS.

Notes:

Message (1) appears automatically.

Error message (2) may appear.

Message (3) appears after a Y to question (1).

B-21
TSPLOT

The TSPLOT program will display contour plots of substrate temperatures on the screen. It is an

interactive program which uses the plot file from SUBSTRAT as input. Because it is interactive and

graphic, it uses certain commands which restrict its operation to IBM PC compatible computers. A plot

file contains the substrate surface (Z = 0 plane) and center-plane (Y = 0 plane, see Figure B-2) temperatures at various times set in the input data file.

Three special files are used by TSPLOT. The CHRSET.VGA file contains the bit patterns for the

graphic display. This file is required. The help file, TSPLOT.HLP, contains the text of the interactive

help messages. Help is activated by pressing the Fl function key. If the help file is not available in the current working directory, no interactive help will be available. The configuration file, TSPLOT.CFG,

sets the colors of the display. The file included on the distribution diskette assumes that a standard VGA

monitor is being used. If the configuration is not in the current working directory, a set of default colors

will be used. These colors may not provide satisfactory contour plots. A new configuration file can be

made by using the MAKECFGG program. See the README file for instructions.

The operation of TSPLOT is briefly explained on the following page. Because the screen interface is similar to TSDATA and relatively few options are available, a detailed description should not be

necessary. After reading this description, try using TSPLOT with the sample plot files. Begin the

program by typing TSPLOT. Abort Uie program by pressing CTRL and C.

Sample Run

At this point the following files should be available in the current working directory: TSPLOT.EXE,

CHRSET.VGA, TSPLOT.HLP, TSPLOT.CFG, TESTl.PLT, and TEST2.PLT. If they are not, use the

INSTALL procedure. If you ran TESTl or TEST2 and aborted before completion, use INSTALL to replace the plot files with the original complete versions.

Type TSPLOT. Press ENTER until the main menu appears; press ENTER again to get the plot file. Enter TESTl.PLT as the name of the plot file. After returning to the main menu, press ENTER to display the next plot. Note that this plot is for t = 1 second, and the pattern is rather small. After

viewing this plot, press any key to return to the main menu. Move (use cursor keys) to the set display

parameters option and press ENTER. Change the Xmin value from "0" to "10", press ENTER, change

Xmax from "60" to "50", press ENTER, and press ESC. Now press ENTER to display the prior plot.

The limits of the display have been changed, and the heated area appears larger. Press ENTER to alternate between the main menu and the next plot. In the later plots note the discontinuity in the profiles

in the padding caused by the air gap. After the last plot you will arrive at the exit option. Instead of

exiting, move to the get plot file option, press ENTER, and enter TEST2.PLT as the name of the plot

file. Return to the main menu and display the next plot which is the first plot from TEST2.PLT. Note

that the limits of the X-axis have not been reset for this new plot. As you continue by displaying

successive plots note how the point of peak temperature is moving along the X-axis. When you reach

the exit option, press ENTER, and respond Y to exit the program.

B-22
TSPLOT Main Menu

This is the main menu screen. It directs you to the different options.

TSPLOT progrcun: ^ Get plot file ^ Display next plot ^ Redisplay ^ Set display parameters ^ Exit

Use cursor keys to move between menu selections.

Press ENTER to activate the menu selection at the X.

Press ESC to return from a selection. Press Fl for help.

The first option is "Get plot file". You must get a (enter the name of an existing) plot file before any

plots can be displayed.

The plot file is read sequentially. You may display the data at the next plot time (option two), or you

may redisplay the last plot shown (option three).

The graphic display consists of three parts: the temperature scale (in °C) to the right, the substrate

surface temperature contours in the upper left of the screen, and the center-plane temperatures in the

lower left. In other words, the temperature contours for the Z=0 and Y=0 planes are displayed together

as if they have been folded along the X-axis so as to both be flat on the screen. The display also shows

the time and the X-axis with coordinates at the left and right edges and tic marks every millimeter.

The fourth option lets you reset the display parameters. You may set the coordinates of the left and right

edges of the display. These values change both the position and the scale of the regions being displayed

thus enlarging or shrinking the plot. The initial limits of the X-axis are for the entire region simulated.

You may also change the temperature contours which are initially set at 25 °C intervals. Changes must

be made so that each temperature is always less than the one above and more than the one below. You

will then probably want to redisplay the last plot.

After displaying one set of plots, you may get a new plot file without exiting this program. This is the

only way to go back in time on the plot file: get the plot file again and start at the beginning.

When you reach the end of the plot file, TSPLOT automatically takes you to the exit option.

B-23

APPENDIX C

ANALYSIS OF NUMERICAL ERRORS PRODUCED

BY A RUNAWAY REACTION RATE

We have made an effort to make a fine grid in the region where the fluxes and flux gradients are high.

Nevertheless, when the reaction rates become very high, the gradients will become so steep that the

approximation of constant temperature and constant reaction rate within a cell becomes questionable. In

this appendix, we examine the magnitude of the errors thus committed by discretization. This analysis can also serve to make appropriate corrections in the program; this has not been done here, partly because

the analysis should first be generalized to the non-symmetric case. (See the assumptions made just below

equation (CI).

Consider the heat diff\ision equation, equation (95). Suppose we have the (correct) temperature

distribution T(x,y,z,t), with a peak at (Xo,yo,Zo); the reaction rate is given by an expression such as equations (108) or (110). Let us simplify this form and assume that

Rp = R„exp(^-TJT(r)) (CI)

where all the preexponential factors are lumped together as the factor R„ and r = (x,y,z) is the position

vector. Consider the terms in equation (12); we simplify the analysis by assuming that the temperature

peak lies at the cell center {i.e., at r^), and that the distribution is symmetric fore-and-aft. Then, since

calculating Rp in that cell means calculating it at the center, and therefore the peak, that means we

overestimate the reaction rate, since the rate falls off" at the faces of the cell.

Symmetry implies that

v(Kvr) = K — + + = 3k

A second-order approximation to this derivative is

V(kV70 = 3k

Because of the assumed symmetry, Ti_i = Tj+j, so that we finally have

T -T

V(kV71 = 6k

Inserting this into equation (12), that equation becomes

{LxY

AT.

pciJ.k Af

= HR(r)- 2k

(Ax)2

(C2)

(C3)

(C4)

(C5)

minus similar terms in Ay and Az (as indicated by the ellipsis). Suppose further, for the sake of

simplicity, that the temperature profile is Gaussian:

e(x,y,x,0 = e,exp

ix-xj'^ + (y-yj'^ + (z-zy (C6)

C-1
where 6 ^ T - T^

and = ambient (reference) temperature.

If we redefine the origin to be at (Xo,yo,Zo), then

Q(x,y,x,t) = e„exp(-r2/02) (C7)

where r^ s + + z^.

Then if x^/a^ < < 1 everywhere within the cell, we can expand (C7). At the center of the cell,

= e(ro) = = (C8)

and

e.-i = e„exp[-(Ax/a)2] » e„[l-(Ax/a)^]. (C9)

Define $ = it^xjof . Then

^.-T-M =e,-6,.i «e^$ (CIO)

so that equation (C5) becomes

(Ax)2 pc— « HR- (Cll)

Assuming three-fold symmetry, the other two terms in (Cll) are the same, and thus

pcM . hR-6k^ (C12)

Af ' " xP- Note that the dependence on Ax has dropped out. Indeed, the last term is the exact diffusive loss rate, resulting from the distribution (C7).

Now let us examine the magnitude of the error made by using a finite value of Ax. From equation (C7),

A0 = 0.-e..i = e„[l-exp(-0] (C13)

Then expanding to the next higher order than was done in (C9),

AO = e„S(l-$/2) (C14)

and the second term in (C5) becomes

^liM . ^(1-5/2) (C15)

(A;c)2

Thus the numerical calculation underestimates the loss rate. Since the correct heat loss rate is

££ = 6x6^/0^ (C16)

and we have

C-2
we multiply L by the ratio

F =

6x6

L m =[l-exp(0].

i . 1 [l-e3q)(-01 1-5/2

(C17)

(CI 8)

in order to get the correct loss rate ^.

As the temperature "runs away," the distribution becomes increasingly peaked, and a declines. As it does

so, ^ increases, and so, therefore, does F; this helps to slow down, or moderate, the runaway. ^ is readily found to be

i = in(e,/e,.i)

Next, consider the source term in equations (C5) and (C12), again. Define

/(W) = exp[-r^/r(je.y,z)]

The correct power output in the elemental volume AV is, from equation (CI),

Lxll Lyil A2/2 R^LV =%R^fdx [

dy //(w)dz

(C21)

(C22)

(C21)

where we have again translated (x^yy^^z^) to be at the origin, and have taken advantage of the 3-fold

reflection symmetry. Again assume equation (C7) to hold. If { < < 1, then within the cell we can

Taylor expand and write

/(x.y.z) =c^(-r^/r.).af(^j^*^

plus similar terms in y and z. Since the peak is at the origin.

and dropping the higher-order terms.

f(x,y,z) « txpi-TJTJ + —

We find that

^{dx'k May'Jo Maz'jo

2T^ 6

(C22)

(C23)

(C24)

/o

(C25)

Thus

f(x,y,z) « txpi-TJTJ

and integration over the cell yields

1 _5i?is.f2 (C26)

C-3
(C27)

This is only valid for Ax/a < < 1. In that limit, however, (C26) is well approximated by

Rp « /e.exp(-r4/r^)exp

which is a reasonable first approximation to the correct integral, even when Ax/a is not very small

The numerical approximation is

4Tj

(C28)

Thus (C28) shows that we must multiply this approximation by exp

better approximation to the source term.

4Tj

(C29)

in order to get a still

C-4
APPENDIX D

ANALYSIS OF IGNITION EXPERIMENT

Among the things we wish to extract from this experiment is the heat transfer coefficient to the surface,

h. Figure D-1 shows the total heating fluxes to ttie gauge, as a function of the power input, with and

without a purging flow of nitrogen. Consider the 18 kW/m^ case. Without the purge flow, the heating

flux is purely radiative, and is measured to be <l>jg^* « 10 kW/m^. With the purge flow on, there is considerable convective heating of the surface as well, and the measured flux is <l>^^i « 18 kW/m^,

suggesting that the convective flux is (^^on ~ ^tot ~ ^rad = 18 - 10 = 8 kW/m^. However, when the

purge flow is turned on, it cools ofl" the heating element somewhat. How much it does so, however, is unknown. Assume that the purge flow reduces the hot-surface temperature such that ^^ad *s reduced to some fraction F of the original. Then

4,^ » lOF (Dl)

and

4),^ = 18 - lOF (D2)

At the point P, on the heater axis but at the surface being heated (5.4 mm below), the heater subtends

a solid angle such that the view factor is Q. Hence the impinging radiative flux is

4)^ = Q€^ar/ + (1-Q)ar/ (D3)

where is the emissivity of the device surface. For the particular case that was carried out, the standofl" distance is 5.4 mm, while the diameter of the glowing filament is 13 mm. According to Siegel and

Howell (1981), Appendix C, the view factor for a disk of radius R, at a point a distance H along the axis

normal to the disk, is

Q =

( W2 (D4)

Hence Q = 0.592. We also assume that = 25 °C « 298 K, and that « 0.85. It follows from

equations (Dl) and (D3) that T4^^. (^-Q)^/

(D5)

If F = 1, this yields T^, « 493 "C, while F = 0.9 yields T^, = 472.5 °C. The device surface cannot

be much colder than 500 °C, since it is observed to glow red. If T^ = 472.5 °C, then (f)^^ = 9

kW/m^.

We will assume that the factor F remains constant for all irradiations. Thus for the case <^tot ~ 25

kW/m^, <l>^^* « 16 kW/m^, and equation (D5) implies that T^ = 840.3 K = 567.1 •'C.

We proceed the same way for all four cases, and obtain the values in the first six colunms of Table D-1

(with fluxes given in kW/m^). Colunm 6, marked 0^, is the convective flux, found as the difl'erence between <^tot and ^^ad- These values of are plotted ys ^tot ^ smooth curve passed through these

points, including the point at the origin. That yields the smoothed convective fluxes given in colunm

7.

D-1
EFFECT OF PURGE FLOW ON HEAT FLUX VS HEATER POWER

MEASURED AT PEAK POSITION; 5.4 MM BELOW HEAT SOURCE

5

4 PURGE FLOW =

0.5 L/MIN

NO

PURGE FLOW

15 20 25 30 35 40

POWER TO HEAT SOURCE (WATTS)

Figure D-1. Total heat flux impinging on gauge, with and without purge flow

D-2
Table D-1. Measured and Inferred Quantities From the Ignition Experiments

Td(K) TdC'C) <*rad hf Tg(K) h f

18 10 745.7 472.5 9.0 9.0 8.9 19.67 698 22.0 0.893

25 16 840.3 567.1 14.4 10.6 10.7 19.56 760 22.9 0.854

34.5 23 921.0 647.8 20.7 13.2 13.0 20.71 837 23.9 0.865

44 32 1000.9 727.7 28.8 15.2 15.2 21.48 907 24.8 0.866

The convective flux can be written in the form

(D6)

where Tg is the purge gas temperature, and T^. is the (cold) gauge temperature. In order to make

progress we make one further plausible assumption: the purging gas takes up a (constant) fraction of the

total energy delivered to the device. We can express this as

where 9 is the temperature referred to the gauge temperature:

Thus, equations (D6) and (D7) yield

(D7)

(D8)

(D9)

Since we have T^ (column 4) and <t>^ (column 7), we readily infer the factor fh from equation (D9). This

is given as column 8.

Next, we must find the heat transfer coefficient h for stagnant flow. The purge air comes down, strikes

the fabric, and must move away radially. This configuration is approximated by the standard problem

of stagnant flow. The heat transfer coefficient for diis case is given in Kaka? et al. (1987). Their

equation (2.176) on page 2-59 gives the Nusselt number:

Nu

= 0.767

\0.43

(DIO)

which was found by Cohen (1961). The subscripts w and e, above, stand for "wall" and "edge" (of the

boundary layer), respectively. The second factor on the right-hand side is approximately one, in this case. The Reynolds number is

Re = (Dll)

D-3
where is the characteristic length. One might think that the characteristic length here is the diameter

of the opening. In fact, however, Cohen (1961) gives Re in the form

^

(D12)

Since the (vertical) velocity goes from to zero in the distance 6, we may write

00 09 dx 8

' (D13)

then equation (D12) indicates that the proper f^. to use here is 6. u^o is readily inferred from the

volumetric flow, dV/dt = 0.5 1/min. Thus

Pa 47;

which yields

= 6.28(y rj cm/sec (D15)

The Nusselt number is given by

iVtt = M (D16)

where 6 is the stand-olf distance between the heater and the substrate, 5.4 nmi. Combining these, and

knowing «(T) and p(T) for air, we obtain h(Tg), as shown in Figure D-2. Note that this is independent

of any estimates of F, the radiative and convective fluxes, etc. We now proceed as follows: we guess a value for Tg; corresponding to this we have h(Tg), and we then

find h(Tg)(Tg - TJ = <^*. We must do this until <j>* = as given in Table D-1. This procedure

converges quite rapidly, and we find the values of Tg and h(Tg) shown in colunms 9 and 10. Finally,

dividing hf by the now-known values of h, we obtain 3ie values of f shown in column 11. There are two

things to be noted about f: it is almost as large as it can get (i.e., close to unity), and it is approximately

independent of </>(oj.

D-4
CNJ CM CM CnI (N CM (N

[>|-3**lu/m] M

Figure D-2. Calculated heat transfer coefficient, as a function of the temperature Tg of the impinging purge

gas jet

D-5

APPENDIX E

THERMOPHYSICAL DATA FOR COTTON

The fabric is a weave of cotton fibers. Cotton itself is a flexible, hollow ttibe of cellulose; the central

channel is called the lumen, and occupies a fraction ri of the total volume. Thus we would expect that the density of cotton is about

p(cotton) = (H)P8,

where Pg is the density of the solid Oargely a-cellulose). Measurements yield t; = 0.2 - 0.4. However,

Brandrup and Immergut (1989) gives the following data on page V-122. (The references given there to the original authors are omitted here for brevity):

Material Density (g/m()

Thus, the density of cotton appears to be very nearly the same as that of the solid (cellulose) implying

that rj ~ 0. We will hereafter ignore the apparently small difference between a-cellulose and cotton, and

take the (mean) density to be p^ = 1.565 ± 0.02 g/mt.

For many polymers, there is only a weak dependence of the thermal conductivity on T, between 100 and

300 K (see Siegel and Howell (1981) Figure 68). One can get an idea of the variability of the thermal

conductivity of cotton and of cotton fabrics from Figures AA-AC, in Kaka? et al.

Measurements at NIST by J.R. Lawson (private communication) have shown that the density of #10 duck

is p « 0.6 g/cm'^, that of #6 duck is 0.72 g/cm^, and that of all the other cotton duck fabrics measured

is

Pf = p(fabric) = 0.62 g/cm^ = 620 kg/m^

The "void fraction" of the fabric is With Pg = 1.565, Pf = 0.62, and with p^ = density of air at STP = 1.774x10"^ g/cm^, the relationship

Cellulose I Cellulose II Cellulose III Cotton

1.582- 1.630

1.583- 1.62

1.61

1.545-1.585

Pf = (1 - 'i»)p3 +

yields $ = 0.6045.

E-1
Next, consider the specific heat. Again, the values given in Brandrup and Immergut (1989) are:

Material Specific Heat (J/g K)

Cellulose

Cotton

1.34

1.22

1.15-1.18 (0-34 °C)

1.32 (0-100 °C)

1.327-1.251

1.273-1.35

Where the temperature range is not listed, it is assumed to be the ambient (20 or 25 °C). We will thus not be far wrong it we take c^ = 1.3 J/g-K,

where the subscript "s" stands for "solid" (cotton or a-cellulose). We expect that, just as was the case

for density,

Cf = c(fabric) = ^c^ir + (1 -

This author has not been able to find the temperature dependence Cg(T); we will defer that question for

the moment.

Finally, we come to the thermal conductivity, k. This is a quantity which is notoriously diflScult to obtain

accurately. Figure E-1, from Childs et al. (1973), shows the large variations in thermal conductivity

depending on measuring conditions. Many of the low values (nos. 5, 6, and 8, for example) were

measured in a vacuum. Numbers 4 and 9, on the other hand, were similar specimens, measured in air

at 25 °C. However, #9 had over three times the density of #4 and would therefore have been expected

to have a substantially higher thermal conductivity.

We list here four sources for the thermal conductivity:

a. The apparently most consistent data from Childs et al. (1973): curves 1, 2, and 3, and

Hence

Cf « 1.122 J/g-K, at T « 300 K.

point 7, give

K « 0.0365 W/m K at T = 306 K = 33 °C

b. We have the following data for cotton paper, from Brandrup and Immergut (1989):

T (°C) K (W/m K)

30

40

50

60

0.049

0.071

0.084

0.090

A curve-fit to these points, if extrapolated, would predict that k vanishes at (and below)

T ~ 10 °C, which is nonsense.

E-2
c. Touloukian et al. (1970) give:

T (°C) K (W/m K)

The density of the material for these

0 0.056 measurements is listed as 81 kg/m^; hence

20 0.058 f> = 0.9493. Then equation (14) gives

100 0.067

K = 0.03408ks + 0.9828Kg (El)

while the tabular data are reasonably well fitted by

k(D = 0.056 + l.lxlO^T, with T in °C (El)

We have K^,i25 °C) = 0.02624 W/m K; now if the data in the table are the conductivities for the light cotton/air mixture, then at 20 °C, with = k^^^, equation (El) yields

Ks = 0.945, which is an unreasonably large value. According to Kunii (1961), the value

for the gas-phase thermal conductivity to be used m equation (14) is 1.7 times what it is in the ambient. If we use k^^ = 1.7x0.02624 = 0.04461, the resulting value of is /Cg = 0.415 W/m K, still very large. If, on the other hand, the listed values are for k^, then = 0.058 and = imply that k = 0.0278, while for = 1.7^^, k = 0.0458.

TTius none of the combinations is plausible.

d. Ohlemiller has made measurements of the thermal conductivity of the (#12) cotton duck;

the apparatus only works properly when the sample is thermally thick, however.

Therefore he carried out a series of measurements, with an increasing number of layers

of the fabric. It was then possible to infer the asymptotic value which would be reached

if the number of layers had been increased to a very large number: It was assumed that

*c(n) = Ko,[l - exp(-ne)],

where 0 = dimensionless thickness of one layer, and n = number of layers. The data were:

n K The Asymptotic value of the series at the left is 0.9.

That is <c« = 0.13 J/m K.

0 0

6 0.050

12 0.077

18 0.099

« = 0.13 is about twice the value (0.056) found in the Handbook, as quoted in c, just

above. Recall that for this fabric, the void fraction is 0.6045. Wiih * = 0.6045,

equation (14) gives

K = 0.28505«8 + 0.84554«

E-3
Figure E-1. Thermal conductivity of cotton as a function of temperature, as measured by different workers

E-4
If we use Jtgas = Kflir, we would then infer that = 0.378, which is an order of

magnitude greater than the earlier estimates. If we use = Uk^^^, then *c = 0.13

implies = 0.324, only a little smaller. The latter value is not too different from the

0.415 found in source c, above.

We thus take = k^(T^ = 25 °C) « 0.324 W/m K

and ie(TJ « 0.13 W/m K

Assuming that the temperature dependence is like that given by eqaution (E2), we find

that

<t(T) = 0.125 + 2.457x10-^ T, (E3)

with T in °C. Alternatively, we might use equation (14) in its more general form, that

is, including the temperature dependence

/c(T) = 0.28505/c,(T) + 0.84554(1.7)/t^i,(T) (E4)

For air, /c(T) is given by equation (52) in Childs et al. For the solid, a reasonable

assumption is that the conductivity is proportional to the absolute temperature, and that

it is also proportional to the density. Thus /Cs(T) is given by

K,(J) = K,(J^ pT/pJ^ s 0.3237 pT/p,T, (E5)

We now turn to the question of the temperature dependence of Cf. For a number of materials, the

thermal diflPusivity is relatively insensitive to T. This is, in particular, the case for wood (see Parker

1985, 1987) which consists, to a significant extent, of cellulose and related compounds. If we assume

that to be the case for cotton, then we could write

/c(T) = apc(T); (E6)

c(T) is generally much easier to measure than /c(T), and this would be a relatively good way to obtain

<c(T); the irony is that we do not have c(T) for cotton or cellulose.

If we take <Cf = 0.13 for the fabric at ambient temperature, we now find that

a s K/pc « 0.13/(620)(1122) = 1.87x10-"^ m^/s

We therefore have, finally,

c(T) = K/pa « /c(T)/620(l. 87x10-'') = 8630*c(T) J/kg K

where /<(T) is given by equation (E3), or by (E4), if preferred.

E-5

APPENDIX F

HEAT TRANSFER COEFFICIENT

IN THE PRESENCE OF A SUBSTRATE

It is first assumed that there is no x-dependence of cigarette or substrate temperature. Now refer to Figure F-1; we need to know how the losses from point P on the cigarette to the substrate surface depend on B. For the convective losses, it is clear that for small 6, they will depend on the substrate temperature

in the neighborhood of P. In fact, for each angle 0, we may identify a mean value y associated with it. For Q = 7r/2, buoyant convection carries gases directly away from the surface, so that the associated

value of y is y = oo , and the convective losses are precisely as they would be for the free cigarette. In

general, we may write

4>,(e) - A,(y)[r,(e) - r/y)] (fi)

where we have explicitly assumed that may vary around the circumference.

It is now necessary to integrate this around the circumference, to obtain the total (convective) energy loss

rate (per unit length). An expression was found relating 6 and y; in order to enable explicit calculations

to be carried out, it was assumed that Tg^y) ^® approximated by the expression

with

Ar= r/y = o) - . (F3)

The integration is carried out in two steps. First, suppose that does not vary around the circumference; the flux can then be integrated, with the result

r\(0)tfe - - Ahie^erfcy/^UT^-TJ

where

\/n ^ Rio, and sfC = RIs

Assuming s = R, a = 0.8R, h^, =10, and Ah = 61, equation (F4) yields

[42.437(r - T) - 20.923 AT]

(F5)

(F6)

Next, the variation in surface temperature around the circumference was taken into account. It was

assumed that

T^iTl/2) - T(l- 20/71)

7',(7C/2)

-7t/2 ^ e ^ 7l/2

Hl2<% <. 3tt/2

(F7)

F-1
where

T = r^(7c/2) - r^(o) (F8)

Note that T^(d=0) is the minimum temperature around the circumference. Hence we may write

7^(0=0) = T,.^in. We will shortly see that t is in fact a function of x and t. The time dependence of

Tg(0) has been suppressed for brevity. Assuming that the used earlier is the maximum temperature,

a simple integration around 6 shows that 7r(18.0366)T/2 must be subtracted from the above result. Hereafter, T^, will always mean the peak temperature around the circumference.

We want to run the cigarette program with model input values of h, T<,, and/or such that the "isolated"

cigarette will have the indicated energy loss via convection. This can be achieved, for example, by using an effective heat transfer coefficient h*, such that

(F9)

0 0

Thus, from equation (F6),

h* - 21.22 - 10.47

( T-T\ s a - 9.02

(Tc -Tc. W/m2K (FIG)
Figure f-i. Schematic of the exchange of energy between point P

on the cigarette surface and the substrate surface.

F-3

APPENDIX G

CALCULATION OF THE

CONVECTIVE VIEW FACTOR

The conductive flux distribution in the transverse direction, </)c(x,y,t), is calculated in Appendix 5-D of

Gann et al. (1988). The i^proximate expression

my,o) « mo,o) (Gl)

is derived there, where 0^(0,0,0) is the (peak) initial conduction/convection flux along the line of contact;

this flux falls rapidly with time. The parameter a is of the order R/2. Equation (Gl) gives ^^.(y)

implicitly.

Experimental measurements (Gann et al. (1988), Section 4) have shown that the (total) flux from a

cigarette can be fitted approximately with a distribution of the form

'^u(x,y,t) = 4>.^ exp

(x-x„-vt \2 _ r (G2)

If expanded as a power series, this agrees with equation (Gl), to lowest order in y^. The total energy

output from the distribution (B2) is

(G3)

Reasonable agreement with experiment is obtained with = 0.6 cm and Oy « 0.32 cm. This suggests

that the above implicit expression for Qc(y), equation (Gl), be modified to be a Gaussian. That is, the

coefficient in equation (139) is

(G4)

where

ho = background « 10 W/m^K and Ah = he - ho « 61 W/m^K.

G-1

APPENDIX H

LONGITUDINAL DEPENDENCE

OF RADIATION FLUX TO THE SUBSTRATE

If the cigarette surface were at the uniform temperature T,., and the view factor for the cigarette for some

point P on the substrate were Q (the fraction of energy emitted by the given substrate area which would

be intercepted by the cylinder), one might think that the radiation flux reaching P is

<|)^ = + (1-Q)<t), (HI)

where

^€.ar/ (H2)

is the flux from the cigarette and

is the ambient flux. If one properly considers the reflections between cigarette and substrate, however,

then with the further assumption that the surface temperature of the substrate, Tg, is also uniform, one

finds that the correct expression would be

_ Q4>. ^(l-OH. ^ a-ec)QH, (H3)

where

,-c,or; (H4)

and

a = (l-£s)(l-£c). (H5)

Equation (H3) reduces to (HI) in the case = I. The more general case, where neither Tg nor are uniform over the surfaces, involves some very complex integrations which can generally not be carried out analytically.

Indeed, is in fact not uniform, but highly peaked. The longitudinal (i.e., parallel to the axis) net radiation flux distribution from the surface of the freely smoldering cigarette,

*,W=..o[r>)-r;] (H6)

can be again be approximated by a Gaussian:

4>^(x) « €,0(Ar)*exp (H7)

H-1
cy

linden

Axis

y

tiai

Figure h-i. The geometric relationship between a general point on the

cigarette and another on the substrate. The profile above the

cigarette is its surface temperature, with its peak at Xq.

H-2
Thus the peak flux from, and peak temperature on, the cigarette surface both occur at x = x^; equation

(H7) gives the peak flux, and equation (H8) defines AT:

r.lx = Tcix,) = + (Ar)* (H8)

Thus, even if we had had the simple approximation (HI) rather than (H3), it must be generalized to

*p = €J

a T^{x)dQ + (1 - Q (H9)

where

jdO, = Q

and where dO is the view factor from a diflferential element on the plane to the difl'erential cylindrical

shell of width dx, a distance x in the axial direction and a distance y in the orthogonal direction, as shown

in Figure H-1. (See reference [20] in Siegel and Howell (1981).) However, we do not have an analytic

expression for dQ.

Two alternative approximations were investigated: in one, the cylindrical shell is approximated (for y >

R) as a wall of height a ~ 2R. In another, the Gaussian distribution (38a) was approximated by a band

of constant (peak) temperature, of equivalent width. That is,

ly^e-^'l'^dx = o^^ (HIO)

This equation yields

a = (Hll)

Thus, if s = ffjj, then with (t,^ = 6.1 mm, a « 1 cm. Neither of these approximations was satisfactory.

Although obtaining the general expression for 0 is a task of daunting complexity, we may obtain a crude

approximation as follows: It is clear that the flux at (x,y) is most strongly determined by the cigarette

temperature at the nearest points. These are just the points at x - that is, on the normal to the cigarette

at X. Carrying this through to its logical conclusion, we may assume that the flux at any point (x,y)

depends only on the local temperature, Tj,(x). This achieves an enormous simplification.

In fact, since each point on the substrate receives radiation from other points on the cigarette as well, the

actual distribution will have a slightly larger half-width than that on the cigarette. The eff'ect can be

expected to be small, however.

Thus at any given distance y from the contact line, the flux distribution over the substrate surface (along

the X coordinate) obtained this way will be similar to the distribution on the cigarette surface, only scaled

down appropriately, as we shall now see. There then remains the task of determining how the flux falls

off" with distance y from the cigarette surface.

From Sparrow and Cess (1978), we find that for a thin but infinitely long strip on the substrate, parallel to the cylinder axis and a distance y from it, the configuration (view) factor Fg^ Q is

= (H12)

H-3
where R is the cylinder radius. This is precisely of the same form as the factor for conduction. (See

equation (Bl).) This is a one-dimensional calculation; clearly, for the two-dimensional case, the flux will

fall off faster. A simple expression which falls off faster but has the same y-dependence for small y is the Gaussian dependence exp(-y^/R^). More generally, the variance of the distribution may be somewhat

different from R (though it may be expected to be proportional to R); call it ffy. Thus we will use

Q(y) = e -yWy (H14)

Since we found above that Oy « 3.2 mm for R = 4mm, we may take, in general,

Oy = 0.8/? (H14)

The improvement (over simply assuming that the heat flux profile is Gaussian in the y direction) which

would be obtained by carrying out the formidable integrations indicated above, hardly justifies the effort. This holds even more strongly for the more complete expression, equation (H3). We thus use

a,(y) = Q(y) = exp(-y2/aj) (H15)

for the sake of simplicity.

Now consider the net flux from the cigarette to the substrate: along the line of contact, 0=1; therefore

equation (H3) yields

^net.rad

I -q ^ '

where a is given by equation (H5).

The radiation flux is treated in the same way as the convective flux, that is,

where is given by equation (142), and this is split into an effective radiation from cigarette to substrate,

,..,

= c,Qo(7-/-r;) (H18)

while the substrate radiates away at the rate

= M *^.a-Q)]o(r/-r/) (H19)

Note that we have taken the reflections from the cigarette partially into account, by using rather than

Cq, however, we have not made the fuller analysis, as in getting equation (H3).

All of the above holds for 0 ^ x < L.

H-4
APPENDIX I

TRANSIENT DEPRESSION OF CIGARETTE

TEMPERATURE UPON DEPOSITION

ON A SUBSTRATE

As noted above, finding T,, njin(t) is a very complex problem; a simple approximation is possible,

however: this problem has been solved for the special case of two semi-infinite solid slabs O^omogeneous, isotropic, uniform) of different temperatures, coming into contact (Carslaw and Jaeger (1959),

Section 2.15). In the vicinity of the contact line, and for a short period, this solution will be

approximately (and locally) valid.

It is assumed there that there is perfect thermal contact between the two surfaces, so that the temperatures

at the interface immediately reach the mean value

T + T

1 +r

where r is the ratio of the thermal inertias of the materials:

r = (12)

N

(*pc).

The temperature of the substrate varies with z and t according to

z

= + Terfc 03)

where z is to be taken as positive into the substrate, and is the thermal diffusivity of the substrate. The temperature of die other "slab" (the bottom of the cigarette) is

7;(z.O = T I + r erf - T.

Since the thermal contact between cigarette and substrate is not perfect, their actual surface temperatures

correspond to small finite depths z, in equations (13) and (14).

For a common fabric/foam substrate, the thermal properties are (Section II.D.2): = 0.1435 W/mK, = 620 kg/m^ and c,, = 1122 J/kg K. Hence = 2.063 x la'' m^/s and the thermal inertia = 316

Ws^^^/m^K. For a cigarette, we find (Mitler, 1988) k^ = 0.316 W/mK, = 620 kg/m^ and c,. =

1043 J/kg K. Hence = 1 • 165 x 10"^ m^/s and the thermal inertia = 293, not very different from that of the substrate. With these two values, we find r = 1.0795. With T^. = 600 °C and T^ = 20 °C,

equation (II) yields <T> « 298 °C.

For this cotton fabric/foam, the temperature fall in a couple of seconds is estimated to be 50-70 "C at the flux peak, at the side of the cigarette; that is. Bid = ir/2. Hence the assumption made above that the

cigarette temperature is unaffected at the side is not quite correct. At 0 = 0, in fact, the temperature fall

I-l
was greater than 200 °C; the precise temperature could not be measured. Thus, ST^^ « 200 °C. Then

from equation (14) we find that erf(x/2^aJ) « 0.38, for t « 2 s. Assuming that the thermal

diffusivities do not vary a great deal among cigarettes, we may take this term to be a constant. In

general, then, assuming that x/l^a^ « const.^ we may write, from equation (14),

6Tmax «r+r-[l+ 0.38 rlf e a * w.^v» • j *

or

« 0.62 rf (15) max

Glass, as a substrate, is a more effective heat sink than a cotton fabric. Indeed, we find that the thermal

inertia for ordinary glass is about 1690; hence in this case, r = 1690/316 = 5.35, and with

T^ = 600 °C and T^ = 20 °C,

and

f = = 97.7 'C

1 +r

« 324 MAX

1-2
APPENDIX J

MASS TRANSFER COEFFICIENTS

Baker and Crellin (1977) measured the diffusion of CO through cigarette paper. For CO in nitrogen, D^, = 0.21 cm^/s, while for oxygen in nitrogen, = 0.199 cm^/s. These are close enough that we may

take the results for CO to be valid for oxygen. They found the relationship:

= AylZ (Jl)

between the diffusion coefficient of CO through the inherently porous cigarette wrapping paper, and the

permeability Z of the paper, to hold. Here:

Dp is in 10'^ cm^/s

and Z is in units of cm^/min (10 cm^ 10 cm water)'^

The value of A is 0.57, in the appropriate units. Permeabilities range between 16 and 2000, in these

units, so Dp ranges between 1 and 22, with many of the papers in the region about:

Dp = 7 X 10"^ cm^/s.

Cigarette paper thickness is 6 = 43 ± 5 /im. Taking the mean, this yields:

7p = Dp/6 =1.6 cm/s

For 7g = the boundary-layer value, we use equation (118). Then at the peak surface temperature

(Tg « 600 °C), the mean value of temperature to be used in equation (118) is 310 °C = 583 K. We

then get:

7j, = 3.3 cm/s.

According to these results, the resistance to diffusion through the paper is only about twice that through

the boundary layer. The combined value, according to equation (122), is

7 « 1.1 cm/s.

M

APPENDIX K

NOMENCLATURE

B = Biot number

Cg = specific heat of the gases

Cg = specific heat of the solid = diflFusion coefficient for oxygen

Djo = dilfusion coefficient for species i in the background "o" [equation (26)]

Dp = mean pore diameter

Ep = activation energy for the pyrolytic reaction

g = acceleration of gravity

h = heat transfer coefficient

hjn = mean convective heat transfer coefficient between cigarette and surface

hq = heat transfer coefficient between cigarette and quiescent air hp = convective heat transfer coefficient; heat of combustion per unit mass

hg = convective heat transfer coefficient (at the end)

h, = heat transfer coefficient for radiation; see equations (15) and (106)

Hj. = molar heat of combustion

k = thermal conductivity of the cigarette (see Section III.D.l, assumption #8)

kg = thermal conductivity of the gas

kop = normalized reaction rate (in min" ^) k, = thermal conductivity of the solid shreds = a characteristic dimension [see equation (34)]

L = length of tobacco column

m = mass of tobacco column

n^^ = no. of grams of ash yielded by the burning of one gram of char

n^ = no. of grams of char produced by the pyrolysis of one gram of tobacco

Hq2 = no. of grams of oxygen which react with one gram of char

Nu = Nusselt number

Pr = Prandtl number

Q,^ = energy released from char oxidation Oower heat of combustion)

Qp = energy absorbed in (endothermic) pyrolysis

Qi = energy loss rates: see Section III.B.7 Q = mean power produced during combustion

r = radial coordinate; ratio of thermal inertias [see equation (12)]; stoichiometric fiiel/oxygen ratio

R = universal gas constant; cigarette radius

Ra = Rayleigh number

Re MqIJv = Reynolds number

Rp, Rgp = pyrolysis rate (in gm/cmh)

Rco = char oxidation rate Sc = Schmidt number (Sc = p/D)

Sj = source/sink term [equation (26)]

T = absolute temperature = cigarette surface temperature at contact line

Tcig = cigarette surface temperature

Tg = surface temperature of the substrate

Uj = radial (convective) gas velocity

K-1
V = total volume of the cigarette

Vg = total volume occupied by the solids in the cigarette

X = coordinate; axis of cylinder representing cigarette

^ox ~ volume fraction of oxygen

y = oxygen mass fraction

Yg = surface value of y

y^ = ambient value of y

y^ = "local" ambient (inside cigarette); used in heterogeneous models

yg = value of y at cigarette (or shred) surface

y-p = tobacco mass fraction

Zp — pre-exponential factor in reaction rate equation (108)

Zqq = pre-exponential factor in reaction rate equation (110)

Greek Symbols

a = thermal dilfusivity of gas

j8 = volumetric coefficient of (gas) expansion; parameter in equation (56)

7 = mass transfer coefficient for air + paper

7b = mass transfer coefficient for air (the boundary layer) alone

6 = characteristic penetration depth

e = emissivity of the fabric (assuming that a, the absorptivity, = c) = emissivity of the (paper) ash

e^ig = emissivity of cigarette

Cj - emissivity of the shreds

Cg = emissivity/absorptivity of the substrate (assumed constant)

6 = e-folding distance in equations (40)-(43)

K = thermal conductivity of the substrate

Kg, Kg = solid and gas-phase thermal conductivities, respectively

V = kinematic viscosity of air Pg = mass density of the gases

p-p = tobacco density

Pg = mass density of the solid a = Stefan-Boltzmann constant T = total duration of the smolder

<l> = void fraction in the cigarette (i.e., the volume fraction of gas, rather than of tobacco shreds) ^ = total void fraction (including the void space in the shreds, and is therefore larger than <^) = ambient radiation flux =

</>c = radiation flux from cigarette; convective flux (substrate)

</>c* = model convective heating flux from cigarette (see equation (9))

0i„ = heating flux reaching the substrate surface

*^net given by equation (6)

<^out ~ surface heat loss of the surface to the ambient, for the areas away from the cigarette

<f>j = radiation flux (substrate)

0 = Q(r) = radiation view factor of the cigarette as seen by the substrate at the point r = fraction giving influence of convection at a given point

K-2
Subscripts

a = ambient; ash

c = convection; cigarette; char

g = gas

m = maximum

0 = original; ambient

p = pyrolysis; peak

r = radial; radiation; effective

s = surface; substrate; solid

t = tobacco; total X = in x-direction (axial)
i
NIST-114

(REV. 9-92) AOMAN4.09

U.S. DEPARTMENT OF COMMERCE

NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY

MANUSCRIPT REVIEW AND APPROVAL

MSTRUCnONS: ATTACH OfUOINAL Of THIS FORM TO ONE (1) COPY OF MANUSCRIPT AND SEND TO: THE SECRETARY, APPROPRUTE EDITORIAL REVIEW BOARD.

TITLE AND SUBTITLE (CITE IN FULt^ Modeling the Ignition of Soft Furnishings by a Cigarette

ERB CONTROL NUMBER

W93-14l7:li

(£RB ySS ONLY)DIVISION

PUBUCATtON REPORT NUMBER

MST/SV 852

PUBLICATION DATl

August 1993

B64

CATEQOflY com

240

mmsBm prsntsd pages

169

CONTRACT OR GRANT NUMBER TYPE OF REPORT AND/OR PERIOD COVERED

Final

AUTHOR(S) (LAST NAME. FIRST INITIAL, SECOND INITIAL)

Mitier, H. E., and Walton, G. N.

PERFORMING ORGANIZATION (CHECK 00 ONE BOX)

NIST/GAITHEnSBURQ

NIST/BOULDER

JILA/BOULDER

LABORATORY AND DIVISION NAMES (FIRST NIST AUTHOR ONLY)

Building and Fire Research Laboratory, Fire Safety Engineering Division

SPONSORINO ORGANIZATION NAME AND COMPLETE ADDRESS (STREET, CITY, STATE, ZIP)

RECOMMENDED FOR NIST PUBUCATtON

JOURNAL OF RESEARCH (NIST JRES)

J. PHYS. A CHEM. REF. DATA (JPCRD)

HANDBOOK (NIST HB) SPECIAL PUBUCATION (NIST SP)

TECHNICAL NOTE (NIST TN)

MONOGRAPH (NIST MN)

NATL. STD. REF. DATA SERIES (NIST NSRDS) FEDERAL INF. PROCESS. STDS. (NIST FIPS) UST OF PUBUCATIONS (NIST LP) NIST INTERAGENCY/INTERNAL REPORT (NISTIR)

i£TTER aRCULAR

BUILDING SaENCE SERIES PRODUCT STANDARDS

OTHER

RECOMMENDED FOR NON-NIST PUBUCATION (CITE FULLY)

| | | | FOREIGN PUBUSHING MEDIUM

PAPER I I CD-ROM

DISKETTE (SPEaFY)

OTHER (SPECIFY)

SUPPLEMENTARY NOTES

ABSTRACT (A 1S00-CHARACTER OR LESS FACTUAL SUMMARY OF MOST SIGNIRCANT INFORMATION. IF DOCUMENT INCLUDES A SIGNIFICANT BIBUOQRAPHY

OR LITERATURE SURVEY, CITE IT HERE. SPELL OUT ACRONYMS ON FIRST REFERENCE.) (CONTINUE ON SEPARATE PAGE, IF NECESSARY.)

This paper describes the user-friendly computer models CIGARET and SUBSTRAT. CIGARET calculates the time-dependent behavior of a cigarette smoldering quietly in the air, away from surfaces. The model

incorporates diffusion and convection of gases, as well as the kinetics of char oxidation. It calculates the

internal heat fluxes, as well as the internal distributions of temperature, gas velocity, and oxygen concentration.

SUBSTRAT determines whether a two-layer soUd (with an air gap in between), exposed to a moving heating

flux such as is produced by a cigarette, will ignite. Among the processes taken into consideration are threedimensional heat conduction in the substrate and its pyrofysis. This model has successfully simulated the

thermal runaway signifying smoldering ignition of the substrate when it is exposed to a set of external heating

fluxes. SUBSTRAT and CIGARET have been designed to work in tandem to simulate the most frequent cause of fatal fires: cigarette ignition of upholstered furniture and bedding. Users' guides are included.

KEY WORDS (MAXIMUM 9 KEY WORDS; 28 CHARACTERS AND SPACES EACH; ALPHABETICAL ORDER; CAPITALIZE ONLY PROPER NAMES)

cigarettes; cigarette model; computer model; free smolder; furniture fires; ignition; mathematical

modeling; modeling; pyrolysis; simulation; smoldering; substrates

AVAiLABIUTY

UNUMITEO FOR OFFICIAL DISTRIBUTION. DO NOT RELEASE TO NTIS. ORDER FROM SUPERINTENDENT OF DOCUMENTS, U.S. OPO, WASHINGTON, D.C. 20402 ORDER FROM NTIS. SPRINGFIELD, VA 22161

NOTE TO AUTHOR(S) IF YOU DO NOT WISH THIS MANUSCRIPT ANNOUNCED BEFORE PUBUCATION,

PLEASE CHECK HERE. | |

ELECTRONIC FORM


IP o

w o

91

3" 03c

TJ W

^'

(0 (fl

8

O 2 C

CD 3 _ 1=1

s a.

ro

- i

w g

{U ft

(DO3-o

o

(a