f‘ 1' 9; 1% I 0 (”Km re‘fliq r». .. g e: -: , 1,5,»? $5., ,.. .M‘RHL - -, -<- A ‘L'x . 9?- . ,- - ,_, g .k .‘ . .j.’ 3'2; «“393 3.15% ‘7; “g, :r 1’ 'W:J.‘r‘-J A“ ' yfl- M vii: 1h. 3.51:”) 3 "T3 . ”'39:“; "5t u Avg‘ :1! {fig (.l inf/W932?“ W '4' -"?;":‘v ’2’“ _ .1. “\ ~ > .‘m(f ’~'.hL4L — _ r 3, b A 'j“"." .-;-x us 31. v.1- L‘ ‘3 Hts-n. “3% flag"; Ii... ‘ fig“. ”,1. .-:uv“ 6;: 1 2:4 if?“ ‘ ‘ : . I, fir: .4233 , . .15. .- ~- - ‘ , . ._.. in: 1.3:,» f“. ;. ti‘l—O " ‘L \n-a, ‘ n We: 9.7:? q - Lax“ - ‘ Huh, “-5": 7 ‘5' ‘- h»m‘?,""‘~."" .mm «m ‘ {'11: ‘ Maggi r .N. w" "'. a? {g}: _ ’fY‘ 3.. m-g 5:}. ,‘l g; .:._. A 3.3—“. urn-3- ' W7. *1}. . m...“ 3£E1§J ”5&1 .¢ u ;£:nu Jr} fi'ermr x 2; 12‘ ’ 5U ‘ . lnnn‘m‘ I~ . 'slmnun DBM‘ I HEM lllllllllllllllllll”(Ulllllllllllllllllllllllllllllllllll 31293 00895 2511 This is to certify that the dissertation entitled THEORETICAL INVESTIGATION OF PILOTED IGNITION OF WOOD presented by LIN-SHYANG .TZENG has been accepted towards fulfillment of the requirements for Ph.D. degree“. Mechanical Engineering ”MC/4a., Major professor Date 3/ 21/] % j/WM 4%%& Co-Advisor [”5 U is an Affirmative Action/Equal Opportunity Institution 0-12771 THEORETICAL INVESTIGATION OF PILOTED IGNITION OF WOOD BY LIN-SHYANG TZENG A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering Michigan State University East Lansing, Michigan 48824 1990 ‘* (6/2311; I 113—- f. (1:, ABSTRACT THEORETICAL INVESTIGATION OF PILOTED IGNITION OF WOOD By Lin-Shyang Tzeng The objective of this work is to understand the details of the piloted ignition process by developing a theoretical model for it. The actual process of piloted ignition of wood is complex. From experimental observations we note that the solid must first chemically decompose to inject fuel gases into the surrounding air, which produce a flammable mixture that may be ignited by the pilot flame. If the fuel evolution rate is too small, only a flash is observed. If the fuel evolution rate is large enough to overcome the heat losses to the surface, then a sustained flame is developed after the flashes. A theoretical model for piloted ignition of a flame in the gas phase above a vaporizing or pyrolyzing solid has been developed. Using this model it has been found that (i) The postulated simplified governing equations adequately explain the pre-ignition flashes that are often observed experimentally; (ii) A rational criterion for .positioning the pilot flame exists; (iii) The heat losses to the surface play an important role, indicating that the fuel evolution rate by itself is insufficient for predicting the onset of piloted i ignition. In this investigation, a numerical integration scheme is. developed that accounts for the often vastly different rates between chemical reaction and convection or diffusion processes in the equations of combustion theory. This new numerical scheme is found to be very efficient for the piloted ignition problem, which involves both pre-mixed and diffusion flames. Finally, a numerical model for piloted ignition of wood which includes transient solid-phase decomposition has been developed. It has been found that the activation energy for the combustion of the evolved fuel is 48 Kcal/mole (also obtained from previous diffusion flame experiments [Puri et. al. 1986]) is more suitable for piloted ignition than 29.1 Kcal/mole (obtained from the curve fitting of the pre-mixed flame data for methane [Coffee et. al.1983]). Also, assuming 100% fuel concentration underestimates the ignition delay, the ignition surface temperature and the mass evolution rate at ignition. However, assuming that the gases evolved from the solid contain 25% fuel and 75% inert produces excellent agreement with the experimental ignition delay, ignition surface temperature and mass evolution rate at ignition. This is also roughly the concentration measured by Abu-Zaid [1988]. Using the analytical solution for wood pyrolysis derived by .Atreya (1983) in the piloted ignition model, the predicted ignition delay is very close to the experimental data obtained by M. Abu-Zaid (1988) for higher heat fluxes (or for heat fluxes larger than 2.5 W/cm2). ii At lower heat fluxes (or heat fluxes less than 2 W/cm2), the piloted ignition model using the analytical solution for wood pyrolysis under— estimates the ignition delay, but the piloted ignition model using the finite difference calculations for the equations of wood pyrolysis predicts the ignition delay very close to the experimental data. This is because the analytical model of Atreys (1983) is valid only in the limit of negligible thermal decomposition, and substantial decomposition occurs during the long time exposures that are characteristic of low heat fluxes. A sensitivity study of the solid phase physical parameters shows that the activation energy, frequency factor and the thermal conductivity of wood are very important parameters for the theoretical model of piloted ignition. The effect of increasing the activation energy is similar to the effect of decreasing the frequency factor; thus, a 10% increase in the activation energy is similar to a factor of 10 decrease in the frequency factor. The effect of thermal conductivity is different from that of the activation energy. For higher activation energy, the ignition delay is larger and the minimum heat flux for ignition is higher. For higher thermal conductivity, the ignition delay is larger, but the minimum heat flux for ignition is lower. iii ACKNOWLEDGMENTS I express my sincerest gratitude to Professors Arvind Atreya and Indrek S. Wichman for suggesting to me the wood ignition area as my dissertation topic and for their encouragement and guidance. Their constant guidance, encouragement were of invaluable help in the completion of the work. I am very appreciate for suggestions and encouragement provided by Professor James V. Beck. I am also very grateful to Professors C. Y. Wang and E. Scott, members of my Ph. D. Guidance Committee, for their helpful comments, suggestions, and encouragements. I would like to thank my colleagues, especially Dr. M. Abu-Zaid and Dr. S. Nurbakhsh for their valuable help and useful discussions. Finally, I would like to thank my wife Wen for typing the manuscript, and my parents for their encouragement. This work was partly supported by the National Institute of Standard and Technology under the Grant No. 60 NANBSOOS78, and my graduate study was partly supported by FIRDI in Taiwan, R. 0. C. iv TABLE OF CONTENTS ABSTRACT ......................................................... 1 ACKNOWLEDGMENTS .................................................. iv LIST OF TABLES .................................................... viii LIST OF FIGURES .................................................. ix NOMENCLATURE ..................................................... xii 1. INTRODUCTION ................................................. l 2. LITERATURE REVIEW ....... - ..................................... S 2.1 Gas Phase Chemical Reaction .................................. 5 2.2 Solid Phase Pyrolysis ........................................ 7 3. A ONE-DIMENSIONAL MODEL OF PILOT IGNITION .................... 10 3.1 Background ................................................... 11 3.2 Mathematical Model of One-Dimensional Piloted Ignition ....... 14 3.3 Numerical Solution ........................................... 18 3.3.1 Premixed Flame Velocity Calculations ................... 20 3.3.2 Numerical Simulation of the Ignition Source ........... 22 3.3.3 Steady-State Diffusion Flame Calculations .............. 23 3.4 Results and Discussion ....................................... 24 3.4.1 Piloted Ignition Phenomena ............................. 25 3.4.2 Quenching of the Premixed Flame (Flashes) .............. 28 3.4.3 Minimum Fuel Rate for Piloted Ignition ................. 34 3.4.4 Location of the Ignition Source ........................ 40 3.5 Significance ................................................. 41 4. EXTINCTION OF A STEADY STATE DIFFUSION FLAME ................. 47 4.1 Numerical Solution for Steady State Diffusion Flame .......... 4.2 The Extinction Damkohler Number for Diffusion Flame .......... 4.2.1 Analysis ............................................... 4.2.2 Calculation of the Extinction Damkohler Number ......... 4.3 Comparison of Results ........................................ 5. THEORETICAL MODEL FOR PILOTED IGNITION OF WOOD ............... 5.1 Background ................................................... 5.2 Mathematical Model for Piloted Ignition of Cellulosic Solids.. 5.3 Analytical Solution for Pyrolysis of Wood .................... 5.4 Numerical Prediction of Piloted Ignition Utilizing Analytical Solid Pyrolysis Solution ..; .................................. 5.4.1 Dilution Study ......................................... 5.4.2 Flashes and Quenching .................................. 5.4.3 Comparison with Experimental Data for Other Woods ...... 5.4.4 Comparison of the Minimum Fuel Evolution Rate with the Fuel Evolution Rate at Extinction of a Steady State Diffusion Flame ........................................ 5.4.5 Effect of Air Velocity and Oxygen Concentration ........ 5.4.6 The Effect of Moisture Content ......................... 5.4.7 Correlation of Results ................................. 5.5 Numerical Solution for Decomposition of Wood ................. 5.6 Sensitivity Study ............................................ 5.6.1 Effect of Solid-Phase Activation Energy ................ 5.6.2 The Effect of Frequency Factor ......................... 5.6.3 The Effect of Thermal Conductivity ..................... 6. CONCLUSIONS .................................................. 6.1 Analytical-Numerical Method for the Chemically Reacting Flow . vi 48 50 55 62 63 65 65 67 75 80 82 88 90 90 95 95 97 99 107 107 110 113a 116 6.2 Analytical Solution for Extinction of the Steady-State. Diffusion Flame .............................................. 6.3 One-dimensional Model of Pilot Ignition ...................... 6.4 Theoretical Model of Pilot Ignition of Wood ............... ... 6.5 Recommendations for Future Work ............................. A. A COMBINED ANALYTICAL-NUMERICAL SOLUTION PROCEDURE FOR CHEMICALLY-REACTING FLOWS .................................... A.1 Introduction ................................................. A.2 The Model Problem ............................................ A.3 The Premixed Laminar Flame ................................... A.4 Benifit of the Analytical-Numerical Method .................. B. PROGRAM LISTING FOR THE PILOT IGNITION MODEL ................. C. PROGRAM LISTING FOR THE STEADY STATE DIFFUSION FLAME ......... D. ESTIMATEION OF SOLID PHASE PARAMETERS FROM PILOTED IGNITION EXPERIMENTAL DATA ............................................ D.l Experimental Data ............................................ D.2 Estimation of Convective Heat Transfer Coefficient ........... D.3 Comparison of the Analytical Solution and Experimental Data... D.3.1 Determination of Thermal Properties of Wood ............ D.3.2 The Effect of h on the Prediction of Thermal Conductivity of Wead ................................... D.3.3 The Pre-Exponental Factor Calculation .................. E. PROGRAM LIST FOR THE PARAMETER ESTIMATION OF WOOD ............. LIST OF REFERENCES ............................................... vii 117 117 118 125 130 I32 189 192 193 195 196 Table Table Table Table Table Table Table Table Table Table Table Table Table Table 5-2 5-3 5-4 5-5 A-l D-l D-2 D-3 LIST OF TABLES Flame velocities calculated by using different grid spacing and time steps. (unites: cm/sec) ............... 44 Steady state diffusion flame location (xf) ............. 45 Location of the ignition source and the fuel concentration at this location. All calculations are for T8 - Tco -298 K ..................................... 46 Effect of the higher (48 Real/mole) and lower (30 Real/mole) activation energies. (Assuming the stoichimetric flame velocity - 38 cm/sec). 74 Comparison of the ignition delay (sec) .................. 87 Comparison of the fuel evolution rate at ignition. (mg/cm sec) ............................................. 87 Comparison of the surface temperature at ignition (00).. 88 Comparison of steady state minimum fuel evolution rate with the numerical simulation of the piloted ignition. (pure methane assumption) .............................. 94 Comparison of steady state minimum fuel flow rate with the numerical simulation of the piloted ignition. ( 25$ methane and 75% inert gas) ....................... 94 Thermal properties of wood as a function of moisture content after parker (1988) ............................ 97 Flame - velocity for Le - 0.5, fl - 10, V - Flame velocity, N - Numer of time step calculated, (CPU time, by Digital Micro Vax II, unit sec) ................................ 131 Thermal conductivity calculate from measured surface temperatue ............................................. 195 Comparison of thermal conductivity by different h. ...... 196 Pre-exponential factor for the pyrolysis, obtained by the best fitting of the experimental data .............. 199 viii Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure .Figure Figure 3-1 3-2 3-3 3-4 3-5 3-6 3-7 3-8 3-9 3-10 3-11 3-12 3-13 3-14 LIST OF FIGURES Piloted ignition of solids -- The physical problem .... Surface temperature-time history during piloted ignition of a red oak sample .......................... Piloted ignition of solids -- The model problem ....... Computational grid. ................................... History of the piloted ignition process from ignition (t - 1.089 s) to development of a steady diffusion flame (t - 3.887 s). Note: the scales for Y and Yf are not shown for clarity. They differ fromothe scale for 0 by a constant factor that is 0.250 for Y and 0.167 for Yf ................................. 8 ........ Heat flux profiles (proportional to temperature gradient) for the propagating pre-mixed flame during piloted ignition ...................................... Temperature profiles for the quenched pre-mixed flame (resulting in a flash) ................................ Fuel Concentration profiles for the quenched pre-mixed flame (resulting in a flash) .......................... Oxygen concentration profiles for the quenched pre- mixed flame (resulting in a flash) .................... Heat flux profiles for the quenched pre-mixed flame (resulting in a flash) ................................ Temperature profiles during flashes caused by the fuel flow rate being below the lean flammability limit ...... Fuel concentration profiles during flashes caused by the fuel flow rate being below the lean flammability limit ................................................. Oxygen concentration profiles during flashes caused by the fuel flow rate being below the lean flammability limit .................................... Heat flux profiles during flashes caused by the fuel flow rate being below the lean flammbility limit ...... ix 12 13 15 19 27 30 31 32 33 35 37 38 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure .Figure Figure 3-15 4-4 5-1 5-2 5-3 5-4 5-5 5-6 5-7 5-8 5-9 5-10 5-11 5-12 5-13 5—14 Minimum distance of the ignition source at surface temperature of 298 and 491 K for various fuel flow rate .................................................. 42 Partial burning ....................................... 56 Diffusion flame ....................................... S7 Near equilibrium, diffusion flame ..................... 59 Match with the same slope ............................. 60 One-dimensional model of the piloted ignition of solids ................................................. 68 A comparison of measured and predicted surface temperature showing the effect of radiant emission. Weight loss shown is due to desorption of moisture. E1. 77 Comparison of regular grid and small grid ............. 81 Temperature history of the gas near the surface (Assume fuel is 100% pure) ............................ 83 Temperature history of the gas near the surface (Assume fuel mass fraction - 0.25) .................... 85 Surface temperature-time history during piloted ignition of a red oak sample .......................... 89 Comparison of predicted surface temperatures and measured surface temperatures at ignition ............. 91 Comparison of predicted ignition delay time and measured ignition delay time .......................... 92 Ignition delay time for different air velocities and oxygen concentrations ................................. 96 Ignition delay time for various moisture contents of wood .................................................. 98 Plots of square root of ignition time vs heat flux for different air velocities and oxygen concentrations .... lOO Plots of square root of ignition time vs heat flux for various moisture contents of wood ..................... 101 Ignition delay time for various moisture contents of wood (solid pyrolysis equations have been solved numerically). ......................................... 104 Density of wood at piloted ignition. .................. 105 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure 5-15 5-16 5-17 5-18 5-19 5-20 5-21 D-l D-2 Plots of square root of ignition time vs heat flux for various moisture contents of wood (solid pyrolysis equations have been solved numerically) ............... 106 Effect of activation energy on the ignition delay ..... 108 Plots of square root of ignition time vs heat flux for various activation energy ............................. 109 Effect of frequency factor on the ignition delay ...... lll Plots of square root of ignition time vs heat flux for various frequency factors. ............................ 112 Effect of thermal conductivity on the ignition delay... 114 Plots of square root of ignition time vs heat flux for various thermal conductivity .......................... 115 Schematic diagram of the experimental facility ........ 190 Comparison of measured and predicted surface temperature ........................................... 194 Comparison of measured and predicted mass flux at the surface ............................................... 198 xi NOMENCLATURE A Pre-exponential factor C Specific heat C Specific heat of unpyrolyzed wood C Specific heat of char d Quenching distance D diffusion coefficient 2 2 Ap C h D Damkohler number, D - -——E——e A -fio E Activation energy F Externally applied heat flux h Boundary layer thickness h Convective heat transfer coefficient Le Lewis number, Le - A/pCpD m Mass flux, m - pv M Nondimensional mass flux, M - meh/A q Heat release Q Nondimensional heat release, Q - q/Cp(TF - Tm) R Ideal gas constant -t3(1-0) .R Nondimensional reaction term, R - DY Y exp[ ] o f 0 l - a (1—0) t Time T Temperature xii v Velocity x Position coordinate Y Mass fraction Creek a0 Nondimensional factor, a0 - l - Tm/TF fl Nondimensional activation energy, 8 - floa 8° Nondimensional activation energy, 50 - E/RTF 6 Laminar flame thickness, 6 ~ A/povon e Emissivity 0 Nondimensional temperature, 0 - (T - Tm)/(TF - T”) A Thermal conductivity Aa Thermal conductivity of unpyrolyzed wood Ac Thermal conductivity of char v Stoichiometric coefficient 6 Nondimensional position coordinate, é - x/h p Density pwd Initial density of wood pwf Final density of wood char 1 Nondimensional time, 7 - At/pCph2 w Reaction rate of fuel(- AszfYOEXP[-E/RT]) wT - qw wo - uw ta Stephan Boltzamann Constant - 1.356 * 10.12 Cal/cmzsecoK xiii Subscripts F Flame f Fuel 1 Ignition source 0 Oxygen 3 Surface a Ambient w Wood xiv CHAPTER 1 INTRODUCTION The burning of wood and other cellulosic materials has long been of central importance to fire research because wood is a primary material for building construction. Therefore, a fundamental study to understand the chemical and physical processes that occur during the combustion of cellulosic materials is a significant step toward achieving effective control and prevention of unwanted fires. Cellulosic materials (such as wood, cotton, paper, etc.), when subjected to heating with sufficient intensity and for a long enough duration, will ignite to yield sustained combustion. Depending upon whether the ignition occurs with or without the aid of an external ignition source, the result is accordingly classified as spontaneous or piloted ignition. The objective of this thesis is to study the piloted ignition of wood theoretically. Although there are some theoretical models for the auto-ignition of pyrolyzing materials (further discussion in chapter 5), there is no theoretical model for piloted ignition of wood. Piloted lignition, like auto-ignition, includes the solid pyrolysis and gas-phase chemical reactions. However, the mechanisms for piloted ignition are different from those of auto-ignition. Actually, the detailed mechanisms 2 are still unclear. From the experimental observations of A. Atreya (1983), the important factors for piloted ignition include: (1) pyrolysis of wood, which produces gaseous-fuels (such as CO, CH4 and higher-molecular-weight hydrocarbons) and inerts (such as C02, H20). (2) diffusion and mixing of the fuel with the ambient air, which produces a flammable gas in the boundary layer. (3) pre-mixed flame ignition by a pilot ignition source. (4) quenching by the surface and (5) establishment of a sustained diffusion flame after the pre-mixed flames has vanished. This thesis develops a theoretical model which attempts to explain and to quantify the above experimental observations. The piloted ignition criterion used in the theoretical model is defined as the development of a sustained diffusion flame after the pre-mixed flash. This ignition criterion is in sharp contrast with ignition criteria used previously, such as: minimum char depth, minimum surface temperature and minimum fuel flow rate. It is important to emphasize that the ignition criterion used in this work is not an arbitrary threshold. Instead, it is a natural consequence of the equation solved in this work. This model would be helpful for understanding the important parameters of fire safety (such as increasing activation energy of wood, decreasing pre- exponental factor of wood and decreasing thermal conductivity of wood). In order to solve the mathematical model which includes both the pre-mixed and diffusion flames, a new numerical scheme has been developed. This method is herein called "a combined analytical and numerical solution method for chemically reacting flow"; the method is 3 described in detail in Appendix A. Usually, the numerical calculation of a pre-mixed flame requires very small time steps and grid spacing, but diffusion flame calculations can use a much larger time step, since the changes are slower. Consequently, the numerical calculations for the transient period spanning both pre-mixed flame and diffusion flame requires is very difficult, often causing over-flow for the numerical calculation. The combined analytical numerical method is comparatively stable during this transient period. In chapter two we present a brief literature review for the development of the numerical scheme for pre-mixed flame propagation, quenching, and pyrolysis and ignition of cellulosic solids. In chapter three a simplified piloted ignition model for the gas phase is developed. Here, irregular grid spacing is used for the numerical calculation. The finite difference expressions are the same as equation 3-100 of Anderson (1984). The computer program has been validated by the numerical calculation of the pre-mixed flame velocity, which is well known for CM‘-air mixtures. Also, the program has been checked by calculating the minimum fuel flow rate for the steady state diffusion flame in two ways, and the diffusion flame location in three ways. In this model, the solid phase is assumed to provide a constant fuel mass flux at a fixed surface temperature. (In the real case, both heat flux and surface temperature vary with time.) Although nthis model is simplified, it reveals the basic structure of the piloted ignition process, which includes a pre-mixed flame in the earliest stage of piloted ignition; this flame may be quenched if the fuel flow 4 . rate is too low (lower than required to establish a steady-state diffusion flame) or if the ignition source is too close to the surface (which is the quenching of the pre-mixed flame). This model also shows that as the pre-mixed flames vanish a diffusion flame appears. This investigation also confirms that piloted ignition is related to the extinction of a diffusion flame. The simplified model of chapter three has been published in Combustion and Flame [Tzeng et.al. (1990)]. In chapter four the extinction of a steady diffusion flame is investigated both numerically and analytically. The analytical solution uses large activation energy asymptotics (AEA) to solve for the, extinction limit of the steady-state diffusion flame. In chapter 5 a complete simplified piloted ignition model is developed which includes the details of the solid-phase pyrolysis. In this model, the effect of fuel concentration on the ignition delay time, ignition fuel flow rate, and ignition surface temperature is investigated. It is found that the assumption of 25% fuel and 75% inert (by mass) gives the most accurate prediction of ignition time, ignition fuel flow rate and ignition surface temperature. A parameter study is also conducted, which shows that the activation energy of pyrolysis, pre-exponential factor and thermal conductivity of wood are very important factors for the theoretical model. Finally, chapter 6 summarizes the overall conclusions of this thesis and suggests various possible avenues of future research. CHAPTER 2 LITERATURE REVIEW As already discussed, the model developed in this thesis must include descriptions of the thermal decomposition of the solid to produce fuel gases, mixing of fuel and air in the boundary layer, pre- mixed flame propagation originating from the ignition source, quenching of this pre-mixed flame by heat losses to the surface, and establishment of a sustained diffusion flame in the boundary layer. To the author's knowledge, there are no piloted ignition models in the literature, Thus, this review basically includes discussions of (1) gas phase chemistry, and (2) solid phase pyrolysis. The piloted ignition model essentially is a combination of the above two phenomena. 2.1 Gas Phase Chemical Reaction This part includes ignition, propagation, and quenching of a pre-mixed flame. The laminar pre-mixed flame speed and flame structure have been extensively studied in the past. Spalding [1956] first introduced a time-dependent numerical scheme for the solution of the ‘equations for conservation of mass, energy and different chemical species. He assumed arbitrary initial profiles and then solved the equations using a standard finite-difference method. Difficulties 6 in the numerical procedure arise because the thickness of the pre-mixed flame is very small (about O.l~0.5 mm, William 1985), thus a grid spacing much smaller than the flame thickness must be employed. Also, the chemical reaction rate in the governing equations is much larger than the diffusion rate, making the profiles very steep and forcing the use of small time steps. Thus, numerical calculations of pre-mixed flame propagation requires both small grid spacing and small time steps. A numerical study of laminar flame quenching was performed by Aly and Hermance (1981); they found that the quenching Peclet number.Pe, increases with decreasing fuel concentration. The theoretical results are in fair agreement with the available experimental data. Spalding (1953) found that when a combustible mixture is ignited by heating a slab of gas, the amount of energy added to the gas must exceed a minimum value for ignition to occur. Numerical integration for slabs of various thickness, initially raised to the adiabatic flame temperature, shows that ignition does not occur for slabs thinner than the pre-mixed flame thickness. For thicker slabs, however, a propagating laminar flame develops. This observation is important for numerical simulation of a pilot flame or an ignition source . 2.2 Solid Phase Pyrolysis Combustible solids can generally be categorised as either (1) those which decompose volumetrically, or (2) those which decompose superficially. Wood belongs to the first group, while some polymers (such as PMMA) belong to the second group. In this study, only wood has been considered. Physically, the properties of wood are highly dependent upon its microscopic structure. Wood consists of cells or fibers, whose diameters vary from 0.02 to 0.5 millimeters. The length of the cell fibers ranges from 0.5 to 1.5 mm. The porosity (ratio of the volume of pores to the volume occupied by the cell walls) of real woods lies somewhere in the range 40-75%. (Kanury 1970). Chemically, wood is composed of three major constituients: (l) cellulose (50% by weight), (2) hemicellulose (25%) and lignin (25%) (Stamm 1964). Pyrolysis of cellulose occurs in two distinct path ways (Shafizadeh 1981). At low temperatures (ZOO-280°C) dehydration of cellulose produces dehydrocellulose and water. This process is slightly endothermic (except in the presence of oxygen) and leads to the formation of char, water, and volatile gases such as CO2, CO, and hydrocarbons. The gases evolved in the dehydrocellulose route are 'primarily noncombustible and the char which remains can oxidize through surface combustion. At higher temperatures (280-34000) endothermic depolymerization of cellulose leads to the formation of tar-like products which are highly condensible and constitute the main gaseous fuel to support a gas-phase flame. Pyrolysis of thick samples of wood involves both the physics of heat and mass transfer and the kinetics of chemical decomposition. Consider a thick slab of wood, initially at room temperature, exposed to a highointensity thermal radiation source. The temperature of the solid gradually increases with time, the surface temperature being the highest. Before the surface layer decomposes, evaporation of moisture occurs and a moisture evaporation zone begins to travel into the solid. At later times, the pyrolysis zone begins to develop and then to propagate slowly into the interior of the solid, leaving behind a thermally insulating layer of char. Regarding the mathematical model of pyrolysis, the formation and growth of a char layer which protects the decomposition zone and the virgin material involves many physical and chemical processes, which makes it very difficult to develop a numerical solution that contains all of them. A few mathematical models of wood pyrolysis have been put forward, starting with the work of Bamford et.a1.(l946). They treated wood as a solid of constant thermal diffusivity. Kung(l972) assumed the thermal properties of the solid varied continuously from their values for virgin wood to their values for char. Kansa et.al. (1977) included a momentum equation for the motion of gases relative to the solid and .obtained good agreement with Lee et.al's (1976) experimental data for low heat fluxes. At high heat fluxes, however, poor atreement was obtained; this was attributed to ignoring the effects of structural ,9 changes (shrinkage and cracking) and to the assumption of a single step-pyrolysis reaction. An analytical solution for surface temperature was derived by Atreya (1983) using the integral method. The solution shows excellent agreement with the experimental results. Base on this surface temperature, an analytical solution for the pyrolysis mass flux of wood has been derived by Atreya and Wichman (1989). This equation is obtained by assuming constant wood density. Thus, the solution is only valid in the early stages of pyrolysis. In this research, both the analytical (integral) solution and the numerical solution of the simplified wood pyrolysis equations have been used to describe the solid phase. CHAPTER 3 A ONE-DIMENSIONAL MODEL OF PILOT IGNITION In this chapter a simplified model of piloted ignition is analyzed. The two-dimensional coupled solid and gas phase problem is simplified by assuming that the mass evolution rate from the combustible solid is a constant , and by employing a plane rather than a point ignition source. With these assumption the model problem requires only a transient one-dimensional analysis of the gas phase. The model equations are solved numerically using the fast scheme discussed in chapter 1 and Appendix A. The pilot flame is modeled as a thin slab of gas that is periodically heated to the adiabatic flame temperature of the stoichiometric mixture. The effects of : (i) the location of the ignition source, (ii) the fuel mass evolution rate from the surface, and (iii) the surface temperature of the solid are investigated. An explanation is produced for the pre-ignition flashes that are observed experimentally. A criterion for positioning of the pilot flame is proposed. The minimum fuel evolution rate, by itself, is found insufficient for predicting the onset of piloted ignition; heat losses to the surface plays an important role. Also, the conditions at .extinction of a steady diffusion flame are found to be practically identical to those for piloted ignition. 10 11 3.1 Background The phenomenon of piloted ignition is rather poorly understood. This is evident from the numerous empirical ignition criteria that have been proposed in the literature. Some of these are: critical surface temperature, critical fuel mass flux, critical char depth and critical mean solid temperature. Of these, critical fuel mass flux at ignition appears physically the most reasonable, but critical surface temperature has proved to be the most useful, since it can be easily related to flame spread. The physical mechanism of piloted ignition is quite complex. The solid must first chemically decompose to inject fuel gases into the surrounding air. This produces a flammable mixture (in the boundary layer), which is ignited by an ignition source. A plausible physical configuration for this process is illustrated in Figure 3-l.In the early stages of solid pyrolysis, the fuel flow rate is very small and the fuel-air mixture in the boundary layer is not combustible. As the fuel evolution rate increasses with time this mixture becomes rich enough to allow a premixed flame to propagate through the boundary layer. This premixed flame consumes nearly all the available fuel and is quenched by heat loss to the surface, unless the fuel evolution rate is large enough to replenish the consumed fuel. Thus, either a flash (a quenched premixed flame) or a sustained diffusion flame -(piloted ignition) is observed. Experimental observations of this phenomenon are shown in Figure 3-2. The momentary rise in surface temperature prior to sustained ignition is due the flashes. 12 .aofipoua Huuaasea one a: meadow no sausages vouoawm «in masons Q 003 \. . N . 23.523 w #3 ._. > mmPHNHz mung w M 3:: onH O; pD~—— - pv(Y - Y - f fs), 6x ayo PD— - PVY . (3-6) 0 6x T - T s - and at x - h, t>0; Yf - 0; Yo - You; T - Tm. (3-7) To simplify these equations, p is assumed constant, implying that m - pv is not a function of x; thus, mass conservation is automatically satisfied. Also, A, Cp and D are assumed constants. The above equations are non-dimensionalized as follows: x At m C h o - —————‘3—- , D - P ram-5°] (3-8) RT CP(Tf - Tm) ao-l-TQ/Tf, . and with the definitions am am 620) + M - 2 , 61 66 86 L {0} - l8 and -B(1 - 0) J R - DYfYoEXP o , l - a (l - 0) the governing equations may be rewritten as L{%f}-{:E’}R- The initial and boundary conditions become: at 1 - O , O < 6 < 1; 0 - 0, Yf - 0, Yo - Y at E - O , r > 0; 9 - 08, an —- - M(Y - Y ) , 82 f fs 8Y0 __. - MYO ’ 36 and at 6 - l , T Z 0; 9 - 0, Yf - 0, Yo - Y (3-9) (3-10) (3-11) (3-12) Equations (3-9) along with the initial and boundary conditions are solved numerically by using a finite difference formulation. 3.3 Numerical Solution A nonuniform grid is used for numerical computations. As shown in Figure 3-4, this grid is closely spaced around the ignition source to l9 .35 Hanoauuuaaaoo Tn. .282 modmmbm “M309 film—EVE gag uo 28, we have yf - O, T'- TL + l - (l + §)z. At 2 - ze, yf - yo - 0 and TL + (a - fi)z - TL + l - (l + ,B)zlz _ ze. Therefore, 1 _ _ 2-79 2 - , giving T - T + . For z - z ~ 0(6), 6 1 4-2? 8 L 1 + 3' e _ _ _ 2 (the reaction zone) we put T - Te - ekof - eklfi1 - e k2fl2 +... _ 2 z - 2 Te B, where c - -—— << 1, We first match to the e T a right of 26, with the same slope for z > ze. Let - 1 is a boundary condition. 8.]; '3“ Now matching the slope to the left of ze, we let 1 2-3 -— - k - l; k, B ° ——— l - -l is the other boundary condition. 6» -00 Again matching the slopes to the right and to the left of ze, from Equation (4-6) we obtain: 59 2‘ x J0 f‘ E 4’ ’0 l ........ 'T'e _ __ yo ' ° . (a - I9) I i v 2e: Yf " 0 TL _ ‘E .F z-O . r029" flow 1. I u z-1 Tu Z - 1 yf-O z - 2 reaction zone e if - 7f e 70-0 _. f - 0 TL 2 FIGURE 4-3 Near equilibrium, diffusion flame. 60 .003. 23. on» and. scans To 80on 2 _ 2 3 a _ d 82 Dfo e 4k1 exp[-Ta/T ] 61 e a?” [1-:o12<1360 3 Defining k1 6 and 6 0 gives 2 d 51 -2 d6 Therefore, 2 d 82 -2 d6 where r dfii d? <15: d? 2(01-E>(01+E>exp[-1. l — (reduced Damkohler number) 6 60 + 651 + ... ._ 2 8 ... 4D§o e exp[-wane] 2 _ 2 (4-13) [1-c0 1 _ _ _ 60(----)(fii-f)(fi1+£)EXP[-(ko€+kifl1)l 60 + £61 (p.-E)axrt~60‘1/31. (p. - 0081 + €)EXP[-6o (01+re>1. <4-14) k0601/3. The boundary conditions are: l as E + w, -1 as E -0 - In Linan's (1974) analysis shows that when 60 is too small, there is no solution for the above differential equation; the extinction 62 Damkohler number 6oe is approximately given as 608 - e[(l-r)-(l-r)2+0.26(l-r)3+0.055(l-r)‘]. (4-15) In other words, when 60 < 608, there is no diffusion flame (or the flame is extinguished). If the activation energy and the pre-exponential factor are fixed (or a particular fuel is chosen for calculations), the reduced Dakahler number 60 and the extinction Damkohler number 608 are then functions of the boundary conditions. In the one-dimensional piloted ignition model which is of concern here, when the fuel flow rate too small, there is no sustained diffusion flame after the flash. This may occur because the reduced Damkohler number is too small. The calculations of the Damkohler number are presented below. 4.2.2 Qa12ulasi2n_2f_fhs_Exfi022120.220kehlsr_nsmheri The fuel used here is methane and the physical properties are the same as given by Coffee (1983). We put: I pA - 3.56 x 10 (l/sec), Cp - 0.323 cal/g K, h - 1.5cm, -5 -8 4 A - 6.2 x 10 KJ/m K sec (900 K of air), pA - 5.61 x 10 g cal/cm sec K, qf - 11355 cal/g of fuel, E - 29.1 K cal/mole, R - 0.001983 'Kcal/mole K, E/R - 14673 K, and u - 4. For every M (M - meh/A) we have 1 Y - Y f,o fs ’ ;§(st + You/V). where Y - l for pure fuel, and £5 63 - Yaw qfo o E -— a - , T* - -————— - 35155 Yf,o (K), Ta - ———, TL - TL/T*’ VY C RT f.0 P * _ _ _ _ Z - 3 Tu - Tu/T* (Tu - 298 K), 8 - TL - Tu’ r - 1 -2 ._, l + a _ _ _ 2 _ _ a - fl pA _ Te T - T + _2 a - 2, D - upAY a, e - -—, 2 - th /A e L 1 + a Cpm f'° Ta P f - l - e-£, Also 60 and 60e are given by Equations (4-13) and (4-15) and we use T - 298 K and Y - l. L fs For a given fuel flow rate m, the reduced Damkohler number and the extinction Damkohler are then easily obtained. The extinction Damkohler number is calculated by employing the following procedure. (1) Assume a higher fuel flow rate, (ii) Calculate the corresponding values of 60 and 608, (iii) Check if 60 < 608; if not, reduce the fuel flow rate and repeat step (ii). This procedure shows that when the fuel flow rate is smaller than -5 2 3.127 x 10 g/cm sec, 6 is smaller than 6 , (6 -0.827, 6 -0.86). 0 0e 0 oe -4.3. Comparisgn 9f Results Here the results of minimum fuel flow rate for piloted ignition are compared with the steady state numerical model and the analytical 64 Damkohler number analysis. The following results are obtained. 1. When the fuel flow rate is smaller than 3.13 x 10"5 g/cm2 sec, there is no sustained diffusion flame after the flashes, This result is obtained by numerically solving the piloted ignition equation (see chapter 3, see final entry in Table 3-3). 2. When the fuel flow rate smaller than 2.88 x 10'5 g/cm2 sec, a steady state diffusion flame is not obtained. This calculation was performed for Equation (3-9) without the a/at term. The boundary conditions are given by Equations (3-11) and (3-12). 3. Both 60 and 60e are functions of the fuel flow rate; for higher fuel flow rates , 60 is greater than the 60¢. When the fuel flow rate is smaller than 3.127 x 10.5 g/cm2 see, the reduced Damkohler number 60 becomes smaller than the extinction Damkohler soe’ indicating extinction. From the above analysis, it seems reasonable to say that for sustained pilot ignition, the Damkohler number for the corresponding boundary condition should be larger than the extinction Damkohler number, which gives a conservative estimate for the minimum fuel flow rate . Chapter 5 IHEQEEIIQAL flQDEL Egg EILQIED IGNITION OF WOOD, 5.1 W Ignition refers to the appearance of a flame in the volatile gas stream evolved from a solid exposed to external heating (usually radiative). Depending upon whether the ignition occurs with or without the aid of an external ignition source, it is accordingly classified as spontaneous (auto-) or piloted (forced). Several ignition criteria have been proposed, such as critical surface temperature at ignition [Simms (1963)], critical mass flux [Bamford et.al. (1946)], critical char depth [Sauer (1956)], critical mean solid temperature [Martin(l965)], etc. Of these, critical fuel mass flux at ignition seems to be physically the most correct since it can be related to flammability limits. However, surface temperature has proved to be the most useful ignition criterion since it can be conveniently related to fire spread [Atreya (1983)]. The ignition of cellulosic materials is a complex process involving both the gas and solid phases. When a cellulosic material (such as wood) is exposed to an intense radiant heat flux, the solid 65 66 chemically decomposes to inject fuel gases into the surrounding air and produce a flammable mixture which can be ignited by an ignition source. Therefore, this process involves flames of both pre-mixed and diffusion type. Kashiwagi (1974) suggests a theoretical model for auto-ignition based on the assumption that ignition occurs when the reaction rate in the boundary layer exceeds a value of about lOJg/cms sec. Gandhi (1986) suggests another model where it is postulated that the auto-ignition point is the gradient reversal in the gas temperature profile at the solid-gas interface. Amos and Fernandez-Fella, (1988) suggest a theoretical model for auto-ignition by considering the absorption of radiation by the fuel vapor as a potential source of ignition. Their model also considers the more practical case of the existence of convective flow over the combusting surface. Their analysis of this evolution process provides information about how a diffusion flame is generated from a combustion reaction which has initially premixed character. A similarly detailed analysis of piloted ignition has not appeared in the literature. Recently, Tzeng et.a1.(1990) have conducted an investigation of piloted ignition. This predominantly gas-phase investigation did not consider the transient solid-phase decomposition. However, it showed the development of a diffusion flame in an initially premixed gas. This analysis is presented in chapter 3. The decomposition products of wood are very complex. Schwenker and Beck(l963) have detected 37 volatile compounds, Goos(l952) has identified 213 compounds as products of wood decomposition. The 67 composition of the pyrolysis products is also a function of the extent of decomposition, the initial density of wood, the temperature at which decomposition occurs, the fraction of char formed, the oxygen concentration in the air, the ambient pressure and the sample moisture content. Abu-Zaid (1988) and Nurbakhsh (1989) investigated the composition of pyrolysis products of wood. They found that the combustible portion is only about 25% (mass fraction) of the total pyrolyzed gas. 5.2 Mathematical Model £9: Eilgted Iggitign of Cellulosic Solids Consider a thick slab of wood placed in an air stream under ambient conditions (as shown in Figure 5-1). The bottom face is considered impervious to both heat and mass transfer and the front face is exposed to a constant radiant heat flux. A part of this incident radiant energy is absorbed at the surface. As the solid surface temperature rises, it radiates a fraction of the energy back to the surroundings. Another fraction is convected to the adjacent air. The rest of the energy is conducted into the solid. Upon further heating, the solid undergoes thermal decomposition. The products of decomposition escape from within the solid through the surface into the gas. These pyrolyzates mix with the heated surrounding air in the boundary layer. Conditions at a distance h above the slab are 'maintained constant by a flowing oxidizer stream. A plane ignition source is placed near the porous plate. This ignition source is periodically ”turned on" to test for piloted ignition. Mathematically 2008 .«o 530% “can 0 E O 0950” a 69 this process is described by the following equations. GOVERNING EQUATIONS: a a e uations The chemical reaction is assumed to be a simple second-order, one- step irreversible Arrhenius reaction, Fuel + v02 0 Products + q (heat) The governing equations are Energy equation: a T a T a A a T ”T p—+pv—-—(——)+—. (s-1) 3 t a x a x Cp 8 x CD Conservation of species Fuel: 8 Yf a Y a a Yf p— + PV— - —- (pr —) - 0. (5-2) a t a x a x a x Oxygen equation: 8 Yb 8 Yb 8 0 Yo p— + pv— - — (72Do ) - we. (5-3) 3 t 8 x a x a x where: 2 . w - Ap Yf'Yo EXPI-E/RT], ”T. q 0, (5-4) (.0 - mo, .1 The initial and boundary conditions are: at t - 0, 0 < x < h: (5-5) Yf - 0, Y - Y , and T - T . O 000 co 70 at t - 0, the radiation source is turned on, resulting in production of the fuel. Thus, the boundary conditions after t - 0 become: At x - 0; a Y (5-6) where pv, st, and T8 are determined by the heat and mass balance of the solid phase. At x - h, t > 0, (5'7) Yf - 0, Y - Y , T - T . 0 000 on o ua 0 Mass balance: a M a p w - ——3 (5-8) a x a t Energy balance: a TV a a Tw PW Cw “" ' “‘ ( Aw ————) (5'9) 71 Decomposition kinetics: a pw g—Z— - ~Aw(pw - pwf) EXP(-Ew/R Tw) (5'10) The boundary conditions and initial conditions are: Tw(x,0) - Tw(o,t) - Ton (5-11) Pw(x00) - PW! Mw(°°.t) - o and the net heat flux into the solid is given by a Tw(0,t) E . . -Awf——;—;-—— - e[ F - :(TS - T0) - 0(Ts - Tco )] (5-12) The gas phase equations are simplified by introducing the Howarth 1 x transformation, 6 - -——J p dx , t' - t. They are nondimensionalized poho 0 2 pAt p0 voC ho T-T0° by defining r - -—-2——2, H— P , 0- , where Tf is the Cppo ho pA Tf-Tco adiabatic flame temperature. Also, 60 - E/RTf, a0 - 1-(Tm/Tf), p-n°*a° we (r -r > o - Ac ,, 2n 23300050 I p f a O p o o 9 13(1 - 0) YOEXPl- ]. Finally, after assuming a Lewis f l - ao(1 - 0) 2 number of unity and p D - constant, the energy equation and species and R - DY conservation become [”H’J L Yf - -l R (5-13) Y6 -v 72 a<~> Ma(-> aQ(-) <5-14) + 2 ar 86 as where L(-) - The boundary conditions are: At 5 - 0, r > 1, 0-9, ‘ s 8Yf/8 5 - M(Yf - st), (5-15) OYb/a f - M Yo, . and at 5 - l, r > 1 0 - 0, ' Yf - 0, (5-16) Y - Y , 0 on - where M is determined by the mass balance at the solid gas interface. The corresponding non-dimensional solid-phase equations are presented in the next section. Ga - the od In Chapter 3, the properties of the gas phase were assumed to be the same as that of a methane-air mixture. The activation energy for the gas phase reaction was assumed to be 29.1 Kcal/mole (as determined by 'Coffee, 1983). While this activation energy yielded good qualitative results, the results did not compare well quantitatively with the ignition data for wood. Recently, Puri and Seshadri (1986) have studied 73 the extinction of a methane air diffusion flame. The activation energy calculated from the experimental data is around 45~47 Kcal/mole. Also, Westbrook and Dryer (1981) found that both higher (48 Kcal/mole) and lower values (30 Kcal/mole) of activation energies can predict the premixed flame velocity quite satisfactorily if the pre-exponential ‘ factor is adjusted. Now, since the flame velocity in a stoichiometric methane air mixture is 38 cm/sec, the pre-exponential factor is adjusted to pA - 9.068x10u(sec-1). The rest of the parameters (pA, Cp’ q, Tf, Too) are the same as before (Chapter 3). (pA - 5.61x10-8 cal g/cm‘K sec, - 2230 K, Tm - 298 K). Cp - 0.323cal/gm K, q - 11355 cal/gm fuel, Tf The effect of higher and lower activation energies on piloted ignition has been investigated. The result is shown in Table 5-1. 74 Table 5-1. Effect of the higher (48 Kcal/mole) and lower (30 Kcal/mole) activation energies. (Assuming the stoichimetric flame velocity is 38 cm/sec) Activation Min. fuel flow Flame temperature energy rate for pilot at min. fuel ignition flow rate _5 48 4.36 x 10 1401 K ,5 30 3.34 x 10 1284 K Extensive studies on extinction of opposed-flow diffusion flames by Ishizuka and Tsuji (1981) have shown that near extinction, the limit flame temperature is remarkably constant for most hydrocarbons and is about 1550 K. Thus, it seems that the higher activation energy for methane is closer to the experimental data for diffusion flames. Since piloted ignition essentially represents the minimum requirement for the existence of a diffusion flame, the higher activation energy is used for the rest of this study. With these parameters, the above partial differential equations and the boundary conditions are solved by a implicit finite-difference method. 75 5.3 a cal o ut o o sis of Wood An analytical solution of the solid-phase temperature field, heat and mass transfer during piloted ignition of cellulosic solids was derived by Atreya (1983). This solution is used here and is outlined below. Consider a semi-infinite solid initially at constant ambient temperature (Tm). At t > 0, the solid is exposed to a constant heat flux, F. Surface oxidation and internal heat transfer between the pyrolysis gases and the solid are ignored. The thermal properties are assumed to be temperature-independent and the net heat required to thermally decompose a unit mass of wood is taken to be zero. Since the heat transfer through an inert solid is by conduction only, the energy balance is described by: 2 30w 6 0w pwcwe-— - Awe-2—, for x > 0, t > 0, at 8x where 0w - Tw - Tm, Tw is the temperature of wood. The solid is assumed to be at the temperature of its surroundings prior to the experiment. Therefore, the initial condition is 0 - 0 for t s 0, x > 0. The boundary condition is obtained from the energy balance at the solid surface and 'is given by 80w -A -- - f (0 , T , t) w8x x—O S s m 2 - c[F - fies - 6((9s + Tm)‘ - Too )1, 76 where e - the emissivity (or absorptivity) of the surface and h - heat transfer coefficient/e, a - Stefan — Boltzmann constant, F - the prescribed incident heat flux, 05 - TS - Tm, where TS is the surface temperature, fs - heat flux going into the semi-infinite solid at the surface, and pw, cw,Aw are respectively the density, heat capacity and thermal conductivity of the solid. Because of the highly non-linear nature of the problem an exact analytical solution cannot likely be found. An approximate analytical solution for the above problem has been obtained by Atreya (1983) by the use of integral methods. Without entering into the details [see Atreya (1983)], the time as a function of surface temperature is obtained as o r (250 +r-m(r+ffi—) a (r+2so ) S S S S t - x ——;2 + 3/211‘. - * , (5-17) 2fs fl (2sos+t+,/F)(r-J79') 8135 2 pwchw _ 2 25 2 2 where K - - -—-2—-, r - - (h + 4aT0° ), s - -—— 0Tco , B - (r -4Fs), 3 e 3 * fs 2 and f - —— - F + r0 + $0 . S S S t In the calculation of the surface temperature, a correction for the energy required to evaporate moisture was made. This quantity was estimated by using the constant weight-loss rate and latent heat of evaporation for water. The result was then substracted from the incident heat flux to obtained the corrected heat flux. This procedure .is mathematically correct. The comparison between the experimental surface temperature and the calculated temperature as shown in Figure 5-2, however, is better than expected. 77 .8..me 2225 3on 2m :39? ”Siege—.2 =< 214.603 6.0" a one .xEE, omnon K £02. 00;" o 60.0"... .oEBx ouvu n "263 vow: «2:305 5630836.. .332 953022 new 026268 5.; 026350 was .638. 628268 5.; 00:02:30 60:03:00 2:0 no woman ozom 3_::c_._Eom a .2 «00:23.3 .302an 2a a... new «a. . I. noon“ use .I .m :8an 2. S 0003 co 3:68.398 "mm "23oz «r ..8 comer «a 02.0 0003 no 3025398 "...m .2239: .0 000933 2 26 w. 550:» «no. 222$ .cofiflEo 0.3.08 .o .023 2.: 052650. 8320062 335m 032020 use 02:32: .0 52500.8 4 1.16.3 ms; . $55on 82.. 84 Bo 86 8» on. o 01 _ . _ _ NI. .0 .2 32 .3 0m «w to. 3.2 .2; cl 8. w |.. On— 01 u f 60 m 88; 28:69 \ N5% :..o . 22.25% 3425.56 \\ 78 The decomposition of wood is assumed to occur according to a first-order global chemical reaction with a known and fixed char density pf; the mass balance of the one-dimensional decomposition process is given by ax at The decomposition kinetics are assumed to be described by 39 W’ — - -Aw(pW - pwf)EXP(- Ew/RTW). 8t and the initial and boundary conditions are pw(X.0) - pw. H§(c,t) - 0. To nondimensionalize the equations, a length scale is defined as L - Adatn/eF, which is obtained as the ratio of thermal conductivity and incident heat flux. Also, a diffusion time is defined as td - 2 L /(Awm/pwc). The other non-dimensional variables are defined as follows: * * * * * T -Tw/Tco' p I'pw/pm, x -xw/L, t -t/td, E -Ew/RT°°, * A - Awtd, * * co * 4 6f - pwf/pw, H - --, M - thd/pr' 2 - 0Tco /F, * * Aw - Aw/Awm' 0 - ow/To. Hence, the nondimensional mass balance and the decomposition kinetics equations are given by * * 8M 8p _*"_? 3x at * and 8p —* - - *(p* - 6*)Exr(-E*/'r*). at 79 Upon integration, Q * * * * * * Ms - A (1 - 8 ) EXP[-E (0 + 1)]dx , 0 * where us is the mass flow rate at the solid surface. During the initial * decomposition process, the density is almost constant, and p is nearly unity. Let AA - -(H* + 42), BB - -(25/3)2 * *2 and O - l + AAO + 880 s s s * * * * For small x , 0 w 08 - x ¢8; therefore, the upper limit of integration * * is replaced by Xmax - 08 /¢s. By substituting this into the above equation, integrating, and keeping only the lowest-order terms of the * asymptotic expansion for large E , one obtains [Atreya and Wichman (1989)] . 'k * 'k * * A (1-6 ) *2 * * EXP[-E (1-1/Ts )1 as - -————;—— rs EXP(-E /Ts ) [1 - *2 . (5-18) o E T S S O a t8 5 The solid phase parameters have been determined by the parameter estimation technique from the experimental data for piloted ignition of Abu-Zaid (1988). The calculations are presented in Appendix D. The .physical parameters used for wood are listed below: _ 2 h (convective heat transfer coefficient) - 2.0 W/m K, 80 Aw (pre-exponential factor for wood) - 6.0x107/sec, Aw (thermal conductivity of wood) - 0.167 W/m K (-Aa), pw (density of the wood) - 0.54 g/cm3 (-pwd), Ew (activation energy for wood) - 31 Real/mole, Cw (heat capacity of wood) - 0.33 cal/g K, 8f (char yield) - 0.25, e (surface emissivity) - 0.75, Tco (ambient temperature) - 298 K. 5.4 iv: a ' :- . o, ' . e- a itio; t_ zi : -na tical olid The numerical solution procedure for the gas phase is the same as before (Chapter 3). The surface temperature of the solid wood is calculated by using equation (5-17), and the fuel evolution rate is calculated by using equation (5-18). It is tactitly assumed that the contribution of any gas-phase exothermic reactions to the rise in the solid surface temperature is negligible during piloted ignition. The uniform grid spacing used in this work is 0.055 mm. A smaller grid spacing was also tested; Figure 5-3 shows the comparison the 0.055 mm grid spacing with the 0.027 mm grid spacing. The heat flux is assumed to be 1.88 W/cm2 and the fuel is assumed to be 100% methane. 'Clearly, the temperature fields near the surface for small and regular grid spacings are very similar. The time for start of flashing for the small grid is 64.7 sec, while for the regular grid it is 61.7 sec. The 81 om 0E0 mfiaomm ..l. 0E0 4.22m .. .. ecu as“ as. as sass to gassed 3 meson Gas “is: om 0* ON 0 n n n u — n n Bafll‘v’HEdWHl 82 time for sustained piloted ignition for the small grid is 71.1 sec, while for the regular grid is 70.7 sec. It is noteworthy that the ignition source has been turned on every 0.5 sec; thus the error is within one time period of turning on the ignition source. It is concluded that the grid space used in this study is appropriate for piloted ignition. Furthermore, the most important aspect of piloted ignition is the development of a sustained diffusion flame after the pre-mixed flame. The smaller grid can only predict a more accurate pre-mixed flame velocity. It will not significantly affect the time for sustained piloted ignition but will dramatically increase the CPU time. Figure 5-4 Shows the temperature history of the gas 0.15mm above the surface. The temperature peaks are due to the flashes and their subsequent quenching. Although it is not immediately obvious from the graphs (due to scaling), it is found that for small heat fluxes, the flashing periods are longer, while for higher heat fluxes, the flashing periods are shorter. 5.4.1 him As mentioned earlier, numerous products are formed during the decomposition of cellulosic materials. For temperatures between 280°C to 500 °C, the primary pyrolysis products are water, 002, C0, -hydrogen, methane, ethane, acetic and formic acids, ethanol, aldehydes, ketones, tar and char. Abu-Zaid (1988) have analyzed the time- integrated product mass and composition of pyrolysis products for 0.1- 83 TEMPERATURE 0.3- M 0.2. 5 E E 0.1- 2.0 FW‘: 9'9“ 1'0 2'0 . ; hemp-a TIME (sec) 'm ('00) 0.3- g «2- E S 1 0.1- 9"“ """" 4’0 ' ' 7"!” “3' TIME (see) FIGURE 54 Temperature history of the gas near the surface. (Assume fuel is 100% pure ) 84 Douglas fir. He found that the total hydrocarbons and carbon monoxide are only about a quarter of the total gas flux from the solid surface. Thus in this dilution study, the fuel concentration of the pyrolysis products is assume to be equal to 25%. Figure 5-5 is the temperature of the gas near the solid surface for the dilution study. The behavior is basically the same as that for the pure fuel case, only the time for flashes and ignition is a little bit larger. Table 5-2 shows a comparison of the ignition delay times for the experimental, pure fuel assumption and dilute fuel assumption. From Table 5-2, it can be seen that the ignition delay for the pure fuel assumption is about 20-30% lower than the experimental value, while the dilute fuel assumption is very close to the experimental value. Table 5-3 is a comparison of the fuel evolution rate for the pure fuel assumption and the dilute fuel assumption. The experimental fuel evolution rate is difficult to obtain from the experimental measurment, owing to the weight loss measurement errors. Roughly, the experimental measurment of the fuel evolution rate at piloted ignition is about 0.2 to 0.25 mg/cmzsec. As can be seen, the fuel evolution rate for pure fuel assumptiom is too low, while the result for the dilute fuel assumption is closer to the experimental data. Table 5-4 shows the comparison of the surface temperature at ignition for experimental, purefuel assumption, and the dilute fuel 'assumption. The surface temperature at ignition for pure fuel assumption is lower than the experimental value. Although the surface temperature for the dilute fuel assumption is also slightly lower than the 85 0.3 0.3- T l 0.2- ‘62 “2‘ 5 o. a.“ E o.“ _ — Havana .. -' fl“ . n 0 . 1'0 ' £0 ' so 9"“ to an i . ME mum %9§2nmnou own: 11 "c 9.3- 0,3. TEMPERATURE {3 ‘ 3 ' 1 . 1—Fmend g__ _ ' - ' . I v—r-uuua 10 h 20 0 IO 20 30 40 i0 20 9.2 - . . —r-1mm 0 2'0 ''''''''''' 40 ‘0 09 ! FIGURE 5-5 Temperature history of the gas near the smface. (Assume fuel mass Emotion = 0.25 -) 86 experimental value, it is closer to the experimental value. Furthermore, the surface temperature for the dilute fuel assumption is within the range of experimental error. From Table 5-2 to Table 5-4, it can be seen that the ignition delay time, the ignition fuel evolution rate and the ignition surface temperature for the dilute fuel assumption (25% methane and 75% inert) is close to the experimental data. Although it is within the experimental errors, the experimental results are consistently under- estimated. Further dilution (20% methane and 80% inert; results listed in the last column of Tables 5-2 to 5-4), brings the calculations closer to the experimental values. This amount of dilution is within the error of the experimental measurements of Abu-Zaid (1988). 87 Table 5-2. Comparison of the ignition delay. External heat flux Experimental Pure fuel Dilute fuel Dilute fuel 2 assumption assumption assumption W/cm (25% CH‘) (20% CH‘) (sec) (sec) (sec) (sec) 1.88 114 70.7 99.0 104.9 2.6 47 33.5 44.5 47.2 2.62 53 33.5 44.3 46.3 3.45 33 18.1 23.6 24.6 3.47 25 18.0 23.8 24.6 Table 5-3. Comparison of the fuel evolution rate at ignition. External heat flux Pure fuel Dilute fuel Dilute fuel 2 assumption assumption assumption W/cm (25% methane) (20% methane) 2 2 2 2 (mg/cm sec) (mg/cm sec) (mg/cm sec) (mg/cm sec) 1.88 0.0364 0.140 0.175 2.6 0.0393 0.144 0.187 2.62 0.0427 0.153 0.186 3.45 0.0418 0.155 0.189 3.47 0.0429 0.170 0.188 (The experimental fuel evolution rate at ignition is 2 about 0.2 to 0.25 mg/cm sec) Table 5-4. Comparison of the surface temperature at ignition. External heat flux Experimental Pure fuel Dilute fuel Dilute fuel 2 assumption assumption assumption W/cm 0 o (25% methane) (20% methane) (C) (C) (C) (C) (C) 1.88 376 318 348 353 2.6 399 330 360 366 2.62 380 332 362 367 3.45 383 340 372 377 3.47 385 341 374 377 5.4.2 Elashes and Quenching From the experimental observations, there are usually some flashes which are subquently quenched prior to sustained piloted ignition. The sharp rises in surface temperature prior to sustained ignition in Figure 5-6 are due to flashes. These flashes are also observed in the numerical simulation. The sharp rises in temperature (prior to ignition) in Figures 5-4 are also due to flashes. The time elapsed between flashes and sustained ignition diminishes as the heat flux is increased. Also, the time between flashes and sustained ignition is larger for the dilute fuel than for the pure fuel, i.e., the time interval over which flashes occur is smaller for the pure-fuel case. 89 .2952. 2.5 .32 s ..c 222:»: 3.23 were: beam: 082:0..32882 Satan 06 EOE 3 mic. omh 00h or» 000 Ohm one one are .211 . _ _ _ _ _ ova aces” cos—:2 .o 22...“ I. .. 1 can . . can so u cos—cu. .0 232383. a can coEeu. .0 as: .5 so 930 238362 . 2 35322 . 2 n :3."— 52a can—a 00 e (3e) aunivaadwal 33'45805 l ouw 1 oz. . «50)5 mm... H xgm v20 0mm no 207520— .. oev 3362 32326 . oov 90 5.4.3 Comparison with Experimental Data for Other Woods. Figure 5-7 shows a comparison of predicted and measured surface temperatures at ignition. The experimental data are obtained from Atreya (1983). The solid line is for the pure fuel assumption, while the dashed lines are for the dilute fuel assumption. From the figure it can be seen that calculations for both of the assumed dilution levels are within the bounds of the experimental measurements. Clearly, the pure fuel assumption under-predicts the ignition temperature. Figure 5-8 shows a comparison of the predicted and experimentally measured ignition delay times. Once again, the pure fuel assumption under-predicts the ignition delay time and the dilute fuel assumption is more accurate. 5.4.4 99mpgrigon of the Minimum Fuel Evolution Rate with the Fuel Evolution Rate at Extinction of a Steady State Diffusion Flame. The minimum fuel evolution rate of a steady diffusion flame (i.e., near extinction) is calculated with the assumption of constant piloted ignition solid surface temperature. The numerical solution is obtained by solving the steady state version of equations (5—13) through (5-16). The analytical solution is obtained by large activation energy asymptotics (Linan 1974). Tables 5-5 and 5-6 show the comparison of the fuel evolution rate obtained from numerical simulation of piloted ignition (transient calculation) with the analytical and numerical calculations of the 91 .sofiuwswfi us mouduwuonaou commune consumes use moHSuouomEou oommudu pouownoun mo consummaoo ~-m th 25>) .3: 6m: 0“”. Fun mwN rMN 0.? coo; 8:8 x . . oom 2:000:02 + H .6an s. womm 200 com o .1. mac: 0 H.8... 33 Nos 825 A ... 99: Nod 8.25 .. m 6.3 0.5a U wows moon womm X H s o 11in. j 11kbiu 1&1111 150 MM Hmw.imW.IM%l. .T ,0 Wm Av wmw 29|¢w u ilLuluQ‘ .. w .. .. ulriWN .. .. .. .. .. .. n m 1% - e w .m x 9 none. (>1) “dwei uomubl 92 .meu mwaoe Coauacwa essences new oEHu heave soauaswa couowuoud mo somHHwQEoo w-m. museum Aommv oE: c222,: com one 002 cm 0 F _ _ . _ _ _ h _ b — _ . _ _ — . _ . . Ooo 2:090:02 x - coo; c0300 + - ..o_aon_ s . .60 Ba . .. 2&0} 0 [0.— 902 some 825 ... - cs: .28 225 .. - E3. menu 0 . e. D .. s ION o o x . 09 e. - [com gum/M xnlj meH 93 minimum fuel evolution rates at extinction of a steady diffusion flame. The diffusion flame is adjacent to a surface whose temperature is the nearly constant surface temperature at piloted ignition for the transient calculation. From these tables, it can be seen that the fuel evolution rate at piloted ignition is quite close to the minimum fuel evolution rate of the steady state diffusion flame (i.e., near extinction). Although the values are slightly larger there is a trend toward increase with increase in heat flux, at least for the values in the first column. 94 Table 5-5. Comparison of steady state minimum fuel evolution rate with the numerical simulation of the piloted ignition. (pure methane assumption) External Fuel flow rate at Min. fuel flow Min. fuel flow heat flux piloted ignition, rate (steady rate (steady- same as in Table state analytical state numerical 5-3 calculation) calculation) 2 2 2 2 W/cm mg/cm sec mg/cm sec mg/cm sec 1.88 0.0364 0.0328 0.0333 2.6 0.0393 0.0325 0.0330 2.62 0.0427 0.0325 0.0329 3.45 0.0418 0.0323 0.0328 3.47 0.0429 0.0322 0.0328 (The surface temperature for the steady state solution is the same as for piloted ignition) Table 5-6. Comparison of steady state minimum fuel flow rate with the numerical simulation of the piloted ignition. ( 25% methane and 75% inert gas) External Fuel flow rate at Min. fuel flow Min. fuel flow heat flux piloted ignition, rate (steady rate (steady- same as in Table state analytical state numerical 5-3 calculation) calculation) 2 2 2 2 W/cm mg/cm sec mg/cm sec mg/cm sec 1.88 0.140 0.128 0.132 2.6 0.144 0.126 0.131 2.62 0.153 0.126 0.131 3.45 0.155 0.125 0.129 3.47 0.170 0.125 0.129 (The surface temperature for the steady state solution is the same .as for piloted ignition) 95 5.4.5 Effect gfi AL; Velocity and Oxygen Concentration The boundary layer thickness is directly related to the ambient air velocity. Because the theoretical model in this study is one- dimensional, changes in the boundary layer thickness are the only method of simulating the changes due to ambient air velocity. Thus, an air velocity of 0.1 m/sec is assumed to correspond to a boundary layer thickness of 1.5 cm and a convective heat transfer coefficient of 2 W/mzx. These values correspond to the experimental measurements. An air velocity of 1.0 m/sec will then correspond to a boundary layer thickness of 0.6 cm and a convective heat transfer coefficient of 2 4 W/m K. Figure 5-9 shows the ignition delay time plotted against incident heat flux for two different air velocities and 02 concentrations. The ignition delay for 15% 02 is larger than that for 21% 02 and the ignition delay time for 1.0 m/sec air velocity is even larger than that for 15% 02. These results agree with the experimental measurements of Abu-Zaid (1988) qualitatively. 5.4.6 W The effect of moisture content manifests itself in changing the physical properties of wood and diluting the products of pyrolysis. Thus, 96 .msowumuucoosoo Cowmxo use mowuwooao> use ucouommHe you cawu aeaov cowuwcwH m-m $0 26 3629, so to 68m Aommv 9:: CODE? ONE cm 0% n n — _ - «Osage SEQ ... NMDOHh Nonim\E 2.03.5 + «ceases. :5 p o; ION (ZWO/M)X”U mew WGPIOUI 97 in this theoretical model, the effect of moisture is studied by changing the thermal conductivity, the heat capacity and the density of wood, The physical properties of wood as a function of moisture content are taken from Parker (1988); these are listed in Table 5-7. The fuel concentration as a function of moisture content is calculated by the following expression: pdry ' pchar YfS(M) - st(drY) * ( )- pM - pchar Here pH is the sample density when moisture is included; hence, pugpdry' Figure 5-10 is the plot of ignition delay as a function of heat flux for various moisture contents; as expected, the higher the moisture content the larger is the ignition delay. Table 5-7. Thermal properties of wood as a function of moisture content after Parker (1988). p (density) C (heat A(thermal fuel cgpacity) conductivity) concentration dry 0.54 1.38 0.167 0.25 11% moisture 0.6 1.657 _0.197 0.218 17% moisture 0.632 1.786 0.214 0.204 27% moisture 0.686 1.974 0.24 0.184 5.4.7 Correlatign of Results In Equation (5-17), if the surface temperature at ignition,(fls), is assumed to be constant, then the following relationship between ignition time and the incident heat flux (F) can be derived, 98 .0003 we mucoucoo ousumaoa msouuw> you came heaoe coauficmH oa-m museum 2:00:00 0:329: .6 Worm Aommv 0E: :23: _ . on: COP on a _ . . . r _ p b h h p n n 23202 NR 8320.: at 83202 s: be n +XOD~ oé (awe/Nomi; mew wepwl 99 _05 F - [Ros t ' + L, (5-19) where 2 pcA K - — _?.. (5'20) 3 e ,o 5 For constant K and L, the relationship between F and t ' is linear, 7208 is the slope of the line, and L is the intercept. Here L may be interpreted as the minimum heat flux for ignition, since it is the heat flux required for an infinite ignition delay time. Figures 5-11 and 5-12 show plots of square root of ignition time vs. heat flux. In both Figures 5-11 and 5-12, the intercept is L - 0.4W/cm2, while the experimental data of Abu-Zaid (1988) show that the intercept is about 1 W/cmz. The reason for this discrepancy is because the analytical solid pyrolysis solution (Equation 5-18) is obtained by assuming a constant density for wood during pyrolysis. For low heat fluxes, a lot of char is formed prior to piloted ignition, and the density near the surface is considerably lower. This means that the analytical solid phase pyrolysis equation over-predicts the mass flux for low heat fluxes, resulting in a smaller ignition delay time. This motivated us to find a complete numerical solution for wood decomposition. 5.1%szme The one-dimensional solid-phase mathematical model for the decomposition process studied here is given by equations (5-8) to (5-12). In these equations Cw and Aw are functions of density. The 100 .mCoHumuusoocoo somaxo use no Hafiooao> was ucouomeu now stm use: m> oEHu coauwcma mo uoou oussum mo muoam HH-n Hummus m o 26 362...; E to some: mmeobéxgu 0.8: E802 9m. Nun whN _ b EN OWN 01 N; . _ _ b «0.x 52>: _ .050 «0.2 P a}. tonne Nos 8350. use (go—038) (g'()—):uouc(ewpr uomufil) 101 . U003 MO mucmuaoo mkfiumHOE WSOHHM> you xSHm ummn m> mafia COHuHGwH mo uoou mumsvm mo muoam NH-n HmDUHh 93:00 83362 $0 5th NNEQEXE Sm: 2520; NAM mN ¢.N ON m; N.— P L _ p b _ _ p — — _ — - L — p _ _ 83%: x8 93%: out 83202 x: be DO¥+ Ide (g’O—OGS) (9°Qé)**(ewn uomubg) 102 density of the solid in the pyrolyzing zone is assumed to change continuously from the initial density of wood; pwd’ to the final density of char, At any instant, a partially pyrolyzed element pwf' of wood in this zone consists of char distributed thrOughout the unpyrolyzed active material. Since zero shrinkage is assumed, all densities are based on the original volume of the wood element; thus p - p + p , where p and p are the densities of the active w a c a c wood material and char respectively. At t - O, p and at w - pwd’ t - o, pw - pwf' The densities pa and pc can be described by the following equations: _ (pw - pwf ) pa pwd _ ’ pwd pwf p - p wd w pc - pwf ( - ). pwd pwf while the thermal conductivity of the solid in the pyrolyzing zone is assumed to be given by where Aa and Ac are the thermal conductivities of unpyrolyzed wood and char, respectively. Also, the specific heat of an element in the pyrolysis zone is assumed to be where, Cpa and Cpc are the specific heats of unpyrolyzed wood and char, 103 respectivity. The differential equations (5-8) to (5-12) were solved by an implicit numerical scheme, and the numerical code was checked with the numerical program written by Atreya [1983]. The physical properties used for the numerical calculation are pwf - 0.135g/cm3, Ac - 0.07 W/m K Cpa - Cw’ Cpc - 0.165 cal/g K. The remainder of the parameters are the same as before. Figure 5-13 shows the ignition delay time calculations using the numerical code for the pyrolysis of wood. A comparison with Figure 5-10 shows that for higher heat fluxes, the ignition delay time . is not much different than the previous result (in which an analytical solution for the wood pyrolysis was used). However, at low heat fluxes, the ignition delay time for the numerical wood pyrolysis calculations is almost four times larger than the analytical wood pyrolysis calculations. This is because a substantial amount of char is formed over an extended exposure at low heat fluxes prior to ignition, as already discussed. Figure 5-14 is the plot of wood density at the time of piloted ignition. From the figure, it can be seen that the density near the surface for the low heat flux is much lower than for the the higher heat flux. .0 s Figure 5-15 is the plot of ( ignition time) ' vs incident heat flux. The slope of the line is different for different moisture 2 contents, and the intercept is about 0.9 W/cm . This agrees with the experimental results of Abu-Zaid (1988). In summary, the one-dimensional piloted ignition model presented 104 .Ahaaao«uo65a om>aom Coon o>m£ «nowadsvo mammaouhn oaaomv boos mo mquuCOo endumaoe «soauu> you mafia huaov coauaGwH mH-m mmoch $629.3, U6 c0328 _ootmEjzv EBcoo @3368 U5 wootm .393 9:5 :02: _ 0mm 00m A 0mm OmVN. om? Om: Om - n n — . 0.56.62 “ANN — — b n — 23m_o: at 83%: a: to O>+§E' 0.? (sz/M)xnu 1qu weplow 105 . and a.“ edge boo: 06.935 .8323 e332 3 e83 mo .33qu 3-... ESE AEEV mOOuCDm U002, 9: got. more—3520 0 .v m N P o L . _ . _ p — n — .. o NEE? u 3: 38: ..I mm o Necks? u we: 60: _..- 106 . Ahaawoauoaac Uo>aom soon 26:. 23.3360 393053 caaomv coo: mo mucoucoo eunumaoa maofizg pom “BAN one; m> 05.3 Coauwcwa. mo uoou ouaacm mo muoam mH-m 5593 $6383 to c0528 _ootoEszv 88:8 83%: .5 68m. $80553... 38: €028. own Nhn wMN RN CNN 0.: N1 wd . . Ll. ? _ _ b . . . p _ p . r p . p p F. b 00.0 ) . m U - w [mod 1. . w u 9 ( w .x. * 13.0 \_/ . .O - 9 ( . . a; 83%: NR + ...smo w 83%: at x - _ 83%: N: o .. ..O S . .. {cc rand . 107 here captures all the essential features of piloted ignition of wood. 5.6 W In the above theoretical study of wood ignition, some of the physical properties were difficult to estimate accurately. Thus, the pertinent question is: how will the errors in the various parameters affect the ignition calculations. In what follows, the influence of solid-phase activation energy, frequency factor and thermal conductivity were investigated. 5.6.1 KW To understand the effect of activation energy on piloted ignition, three different values were used in the calculations. Figure 5-16 shows the ignition delay time calculations for Ew-28 Kcal/gmole, 31 Kcal/gmole and 34 Kcal/gmole. It can be seen that a 10% change in the activation energy significantly changes the ignition delay time. -For higher activation energy, the ignition delay time is much longer than that for the lower activation energy. This high sensitivity to activation energy also suggests that piloted ignition measurements may be a very accurate way of determining an overall activation energy for the one-step solid phase decomposition model. .0 5 Figure 5-17 is a plot of (ignition time) ° vs heat flux. It 108 com £316 Coauficwfl as...“ no mwuoco aoauuZuom mo uoommm waé 559B 355 cozo>zo< do womtm Aommv mE: coEc _ 00¢ own. omVN om: . o_oE W\_oox¢n o_oE m\_oov:n o_oE ”<60wa _ - ¥ + D o o; —(zw:>/M)>auou msouuu> w“ mo uoou mumsvm mo macaw 5H-n mmowum 208 "<60an 29: m\_oov_ F, 29: @58me b . + Wm 36: cozo>zo< do bootm SEQ 5x5 38: 35%; mam N6 mN TN o.N m; N; md _ll- _ b _ p — n _ — — - n _ — . . p _ _ _ h — _ OOoO Imod :26 Lewd (9°O—398)(g'O—)**(8w!l uomufiI) 110 can be seen that the intercept is 0.6 W/cm2 for Ew - 28Kca1/gmole; For Ew - 31 Kcal/gmole the intercept is 0.9 W/cmz, and for Ewe 34 Kcal/g mole, the intercept is 1.4 W/cmz. The slope of the three lines is therefore also changed. This is because the surface temperature at ignition (which affects both the heat loss and the slope) is changed with change in activation energy. The surface temperature for the three cases was 350, 410, and 470°C. Since the intercept of the line is interpreted as the minimum heat flux for ignition, the higher activation energy yields a higher minimum heat flux; a 10% increase in the activation energy results in a 50% increase in the minimum ignition heat flux. 5.6.2 WWW In the literature [Atreya (1983)], the frequency factor for the pyrolysis of wood ranges from 6.0x101 to 7.5x108. In this study, three frequency factors Aw-6.0x101, 0.6x101 and 60x107 were investigated. Figure 5-18 shows the ignition delay for various frequency factors. As expected, higher frequency factors yield lower ignition delay times. Figure 5-19 is the plot of (ignition time)-o'5 vs heat flux. The trend and the slope of the separate lines is similar to the activation energy case. The intercept for different frequency factors is different. -For Aw - 0.6x101/sec, the minimum heat flux for ignition is about 1.4 W/cmz, for Aw - 6.0x107/sec the minimum heat flux is 0.9 W/cmz, 7 2 while the minimum heat flux for Aw - 60x10 /sec is 0.6 W/cm . In 111 53% do??? one so uouoam hoaosvoum mo uoommm «Tn 356: combo“... :ocmnomi .6 womtm Aoomv o8: coEc _ own owe own _ cam ems ‘o oom\hm.oo x oom\hmo.m + r 8m\Bo.o e. . o; (aw/mm”. wen wepml 112 .muouoem moaosuoum msoaue> you stM was: m> mafia Goauflawa mo uoou enemas mo muoam ma-m b00000 :ocmsomcu U6 sootm ANEQEXE 38: 3.5%; mwN ¢.N ON 02 N; pr—Lnbr—L— — _ — b _ b — . — o8\B.om omm\nmm aimed D + X gunk 00.0 INm.0 (9°o—OeS)(.9°O—)**(9w!1 U0Hlufil) 113 a other words, a factor of ten reduction in the frequency factor increases the minimum heat flux at ignition by 50%, whereas a 10% increase in the activation energy has a similar effect. This is caused by the linear dependence of the the reaction rate on Aw, and the exponential dependence of the reaction rate on Ew. 5.6.3 113W The literature values for the thermal conductivity of wood range from 0.134 W/mK to 0.263.W/mK [Atreya (1983)]. In this study, the thermal conductivity was assigned the values Aa - 0.25 W/mK, 0.167 W/mK and 0.111 W/mK. Figure 5-20 shows the ignition delay time for these various cases. Although the ignition delay time for lower heat fluxes is very large, the three curves tends to practically the same asymptotic value. Figure 5-21 is the plot of (ignition time).o'6 vs heat flux. For Aa-O.25 W/mK, the intercept is about 0.8 W/cmz, for Aa-0.l67 W/mK, the intercept is about 0.9 W/cm2 and for Aa - 0.111 W/mK, the intercept is about 1.0 W/cma. The thermal conductivity is related to the slope of the line. From equation (5-19), for higher thermal conductivity the slope is smaller and the intercept [or L in equation (5-19)] will be lower. This is confirmed by Figure 5-21. Clearly, the higher the thermal conductivity, the larger is the ignition delay time and the smaller is the minimum heat flux for ignition. Physically, this .means that as the thermal conductivity increases, more heat have been conducted into the deeper area, the surface temperature and the pyrolysis product decreases. 1131) From Figures 5-17 and 5-21, it can be seen that for higher activation energy the ignition delay is longer, and the intercept of the line is also larger, whereas, for higher thermal conductivity, the ignition delay is larger but the intercept of the line is lower. This is because the ignition temperature for higher thermal conductivity is lower (for Ail-0.25 W/mK, 'r ignition is about 400°C, for Ail-0.167 W/mK, T ignition is about 410°C, and for Aa-0.111 W/mK, T ignition is about 420°C). However, the ignition temperature for higher activation energy is larger. (for Ew-34 Kcal/gmole T ignition is about 470°C, for Ew-31 Kcal/gmole T ignition is about 410°C, for Ew-28 Kcal/gmole, T ignition is about 350°C). These changes in the ignition temperature in conjunction with Equation 5-19 explain the observed trends of ignition delay for different activation energies and different thermal conductivities. 114 33503300 _oEtoE. 00 5th Aoomv oEz coEc0_ 39.26 383.2%.“ 05 so huffiuosvaoo Hoe—no.5 mo uoommm 0N-m 959nm one can cam one .o x $263 n .28 x x 55583 u .28 + - x E\:,:.o u .28 o . - 0.— (zwd/M)>«uo.«6aoo Hoe—Hosp mags; pom XSC use: Er 2.3 :03ch.“ mo uoou ouasvm me 30$ 3%. HgGHh 333828 6:88: to 68m. Ameobéxo: so: 38%; 0.». Wm m.m _ .2 9m m; N; we _ b n — n n — p b - PL.- - — p P . — p n — b 00.0 .86 [2.0 Lewd o. {:86 u .28 x . v. Ebébco u .28 + . x E\::.o n .350 o - {no (g'Q—OSS)(9°Q—)**(8LUD, uomufil) Chapter 6 Conclusions This chapter presents a summary of the important results for different components of this thesis. Section 6.1 summarizes the conclusions of the analytical-numerical method for the chemical reacting flow, section 6.2 summarizes the conclusions of the analytical solution for the extinction of the steady-state diffusion flame, section 6.3 summarizes the conclusions of the one-dimentional model of piloted ignition, and section 6.4 summarizes the conclusions of the numerical simulation of piloted ignition of wood. Some recommedations for future work are presented in section 6.5. 6.1. a ca -N erical ethod for the Chemicall eactin Flow 6.1.1. The required CPU time for the analytical-numerical method is less than for the implicit scheme. The explicit scheme for the chemical reaction scheme is fastest, but it is less stable when larger time steps are used. 6.1.2. The analytical-numerical scheme is very efficient for the piloted ignition problem. Piloted ignition involves both pre-mixed and diffusion flames. The pre-mixed flame requires very small time steps, while the time step for 116 117 diffusion flame can be 0(100) times larger than for the pre- mixed flame. The implicit and explicit methods overflow very easily during the transition from pre-mixed flame to the diffusion flame, while the analytical-numerical method does not overflow. 6.2 Analytical Solution for the Extinction of the Steady-State Diffusion Flame, The extinction limit for a steady diffusion flame can be derived by large activation energy asymptotics (AEA), which provides expressions for the ignition and extinction Damkoher number. By using a higher mass flux, the flame temperature and the flame Damkohler number are obtained; then by reducing the fuel mass evolution rate, the Damkohler number is again calculated. If the Damkolher number is smaller than the extinction Damkolher number, no solution for the diffusion flame exists, and the flame will be extinguished. The minimum fuel evolution rate for piloted ignition is very close to the analytical solution for the steady diffusion flame extinction limit (110%). 6.3. One—dimensional Model of Pilot Ignition, 6.3.1 The simple one-dimensional numerical model of piloted ignition describes the basic strunture of the piloted ignition process, which includes a pre-mixed flame, the quenching of a flash, 118 and the development of a diffusion flame after the flash. 6.3.2 The minimum fuel evolution rate in itself is not sufficient to predict the onset of sustained piloted ignition, and the heat losses to the solid surface play an important role. 6.3.3 The optimum location of the ignition source is found to be the eventual location of the steady diffusion flame. 6.4 MW 6.4.1. 6.4.2. Although the thermal conductivity of wood is a function of the exposed heat flux, temperature, moisture content and density, the overall thermal conductivity can be determined by the least square fitting of the experimental surface temperature. The curve fitting uses the analytical solution for the surface temperature of wood derived by Atreya (1983); the calculated thermal conductivity of dry Douglas fir is within the range of values listed in the literature. The integral analytical solution for the fuel mass production rate in the initial stages of decomposition can be used to determine the pre-exponential factor of the pyrolysis of wood by using the least square fitting of experimental measured mass flux around piloted ignition. The predicted pre-exponential factor for dry Douglas fir is within the 6.4.3. 6.4.4. 6.4.5. 6.4.6. 119 range of values listed in the literature. Since the analytical solution is only valid for the initial stage of decomposition, this model cannot be used for lower heat flux where ignition delay is very large and a lot of char has been formed before piloted ignition. By using the thermal conductivity obtained by the least square fitting, the analytical solution can predict the experimental surface temperature of wood very accurately. This numerical model for the piloted ignition of wood can predict the flashes, the quenching and sustained piloted ignition. For lower heat flux, the period of flashes and quenching is longer, while for higher heat fluxes, the flashes and the quenching are not so obvious. During the flash-quenching period, the fuel evolution rate is much lower than the minimum fuel evolution rate of the steady state diffusion flame, whereas at sustained piloted ignition the fuel evolution rate is very close to the minimum fuel evolution rate. This leads us to believe that conditions for piloted ignition for wood are quite similar to those for extinction of a diffusion flame. The assumption that the fuel evolute from the solid wood is the same as pure methane underestimates the ignition delay, the ignition surface temperature, and the ignition 6.4.7 6.4.8 120 mass evolution rate. The predicted ignition delay is about 30% lower than the experimental value, the predicted ignition surface temperature is about 15% lower than the experimental ignition surface temperature, while the predicted mass evolution rate only about 1/5 to 1/7 of the experimentally observed value. The assumption that the fuel evoluted from the solid wood is 25% methane and 75% inert gas (by mass) enables the experimental ignition delay, the ignition surface temperature and the ignition mass evolution rate to be accurately predicted. From the experimental observations, the gas flow out of the solid in the early stages of pyrolysis is primarily water vapor, C02 CO, CH‘, and higher hydrocarbons; it is evident that the concentration of total hydrocarbons and C0 is very low at piloted ignition. Both the surface temperature and fuel evolution rate are very important for the piloted ignition, but the dilusion of the fuel is also a very important factor. The study of the effect of 02 concentration and ambient air velocity shows that higher air velocities and lower 02 concentrations increases the ignition delay and increase the minimum heat flux for piloted ignition. This result is agreement with the experimental observation of Abu-zaid (1988). 6.4.9 6.4.10 6.4.11 121 The study of the effect of moisture content shows that a higher moisture content increases the ignition delay without altering the minimum heat flux for ignition. The model which use the analytical solution for the wood pyrolysis equations [derived by Atreya and Wichman (1989)] predicts an ignition delay at high heat flux (greater than 2.6 W/cmz) very close to that calculated from a model which solves the pyrolysis equations numerically. At lower heat fluxes (less than 2.0 W/cmz), the ignition delay predicted by the analytical solid solution model is only 1/3 of the value predicted by the numerical solid-phase model. The plot of (t ignitionffl'l5 vs heat flux shows that the numerical solid phase model obtains closer value of minimum heat flux for piloted ignition to the experimental observation of Abu-zaid (1988). The analytical solid phase model underestimates the minimum heat flux for piloted ignition. A sensitivity analysis shows that: A. The activition energy for thermal decompostition of wood is very important for the prediction of piloted ignition. An 0(10%) change in the activation energy makes an 0(100%) change in the ignition delay. The plot of (t ignition).o'£5 vs heat flux shows that the curves for different activation energies have similar slopes, but that the minimum heat fluxes for piloted ignition are different. An 0(10%) 122 increase in the activation energy makes an O(50%) increase in the minimum heat flux. B. The frequency factor is very important for the prediction of piloted ignition, but is less sensitive than the activation energy. An increase in the frequency factor reduces the ignition delay. The plot of (t ignition).o 5 vs heat flux shows that the curves for different frequency factors have the same slope, but the minimum heat flux for pilot ignition is different. (This trend is similar to the effect of activation energy). The effect of a factor of 10 increase in the frequency factor is similar to reducing the activation energy by 10%. C. The thermal conductivity is very important for the prediction of piloted ignition; the higher thermal conductivity, the larger the ignition delay. The effect of thermal conductivity on the ignition delay is different from the activation energy. The plot of (t ignition)-o 5 vs heat flux shows that the lines with higher thermal conductivity have a lower slope, and that a higher thermal conductivity produces lower minimum heat fluxes for piloted ignition. 6.4.12 A new definition of piloted ignition is proposed in this work. It is found that the development of a sustained diffusion flame after the pre-mixed flash successfully predicts the experimentally measured ignition delay, ignition surface temperature, ignition mass flux and minimum fuel flow 123 rate for piloted ignition. This model is helpful in clearifying important parameters for piloted ignition that are not easy observed from the experimental data. (such as quenching effect, location of ignition source, thermal conductivity of wood and activation energy of wood). 6.5. Recommendation or the t e Work In the present model of piloted ignition, the chemical reaction in the gas phase is assumed to be a one-step irreversible chemical reaction. In some pratical experimental observations, the one-step chemical reaction is too simple and multistep chemical reactions may produce more accurate results. Besides, the composition pyrolysis products are very complex, a great deal of C0 and water vapor are contained in the pyrolyzate; thus, a multistep chemical scheme would give a better representation for the effect of water vapor and CO during piloted ignition. For the time being, the piloted ignition model in the gas phase is one-dimensional. The actural gas phase in the boundary layer is two dimensional. Since piloted ignition involves pre-mixed flames, and since the pre-mixed flame thickness is very small (only about 0.1mm), the grid spacing should be smaller than the flame thickness. A two-dimensional .model for the gas phase would be more practical, but such a model might require a more powerful tool, such as a super computer. Appendix A A COMBINED ANALYTICAL-NUMERICAL SOLUTION PROCEDURE FOR CHEMICALLY-REACTING FLOWS A.1. INTRODUCTION The possibility of exploiting the often vastly different time scales for chemical reaction [0(10-ssec)] and for the subsequent convection and diffusion of this chemical release into the flow field [0(1 sec)] appears to have been first discussed physically by Baum and Corley(1981). (A more detailed history of the ”operator-splitting" method is given in the recent study of Wichman (1990)). In their short work Baum and Corley (1981) analyzed the laminar mixing of fuel and oxidizer streams in a plane channel flow for arbitrary Lewis numbers, and for £39 reaction models; (1) delta-function heat release, and (2) one-step Arrhenius kinetics. For (2), the solution is obtained by dividing each streamwise marching step into two parts. In the first part, a locally homogeneous rate chemistry problem is solved by ~quadratures (for (l) the analytical solution of this part is trivial). This solution is then applied as a (local) initial condition for the second part, in which this heat release is transported by convection and diffusion to the neighboring fluid elements. Thus, the solution is 124 125 obtained by an "operator-splitting" procedure. The purpose of the present study is to develop a rational motivation for the use of the operator-splitting procedure in Combustion theory by analyzing a simple model problem in some detail. A numerical scheme that is consistent with this model is then applied to a simple premixed laminar flame model that has been previously employed as a test for various numerical integration schemes [Peters (1982)], such as the method of lines, adaptive grids explicit and implicit schemes. A.2. The node; Problem Consider a purely thermal model, for which the gas-phase energy balance is represented by transient, diffusion and heat-release terms. When properly nondimensionalized, the coefficients of the transient and diffusive terms are unity, while the heat-release term contains the factor D - Atd/(Atc), called the Damkohler number. Here Atd is the characteristic diffusion (or convection) time, while Atc is the characteristic reaction time. Thus, we consider the following initial boundary-value problem, 2 61' 31: ] —--—2+DQ(x,t), -oo0 at 6x l (M) T(“D,t) - T(Q,t) - 0’ T(x,to) - To(x). J Note that the heat-release term (analogous to the chemical reaction term 126 in the combustion equations) is written as a simple function of x and t. The solution for T(x,At) can therefore be obtained from the solution at the previous time step to from the equation (we put to = 0 with no loss of generality) m At“) T1dedr - J To c - - EXP[-————————-1 (A-3) 2J«(At - r) 4(At - r) is the Green's function for equation (A-l). The first term in equation(A-2) describes the influence of the inhomogeneous reaction term, while the second term describes the redistribution of energy from the previous time step by the homogeneous diffusion operator 2 2 8(0)/6t - a (0)/ax - 0. Note that the second term does not contain Q. The difficulty with obtaining numerical solutions of equations such as equation (A-l) arises from the largeness of D. This usually necessitates taking very small time steps, of order Atc. Rational means for avoiding this requirement (and therefore speeding up the computation substantially) are sought by analyzing equation (A-2) as D 4 w. In order that all three terms in equation (A-2) balance as D 4 m, the time interval At in the first integral must become very small; in _1 fact, At must be of order D z Ate. Thus, by defining s = Dr as a new 0(1) variable of integration, one finds, as D 4 w, 127 CD C(X.AtI£.0)Q(€.0)d€ - J Toc> 0(1), as required for the anaylsis of Sec. II. A.4. Benefit of the Anglytical-Numerical Method The numerical integration scheme of Sec. III was compared to the flame-speed calculations presented in the GAMM workshop [Peters (1982)], and to calculations of the same problem using an explicit numerical scheme1 and an implicit numerical scheme (iterate with the reaction term). An example of these latter comparisons is shown in Table l for the case Le - 0.5, fl - 10. The numerical scheme of Sec. III converges to a solution in all cases, whereas the other schemes develop convergence and overflow problems. The required CPU time is less than for the implicit scheme, though the explicit method is fastest when it works. In multiple space dimensions it is believed that this method will be much quicker and easier to implement than other, more sophisticated algorithms. In addition, flames of mixed (diffusion/premixed) character are not expected to produce excessive computational difficulties. Remarks: 1. The scheme of Sec. III is essentially "adiabatic", while the explicit scheme is "constant temperature". Therefore, they represent, in a sense, thermodynamic extremes. Explicit schemes have also been used by Reitz(l981). 131 Table A-1: Flame- Velocity for Le-0.5, fi-lO, V-flame velocity, N-number of time step calculated, (CPU time, by Digital micro VAX II, unit sec), t - 0.001 0.002 0.005 0.01 0.02 0.05 ‘ 0.1 method 1, case 1 V-0.962 V-0.960 V-0.945 V-0.944 V-0.927 V=0.850 V-0.624 N-900 N-450 N-150 N-lOO N-120 N-40 N-40 CPU 173.04 88.54 32.17 22.98 46.28 37.79 70.54 Ax-0.25 method 2, case 1 V-0.965 V-0.965 V-0.959 not N-900 N-450 N-150 convergent CPU 204.66 172.91 162.76 method 3, case 1 V-0.964 V-0.963 V-0.954 V-0.938 V-0.952 over N-900 N-450 N-150 N-240 N—SO CPU 39.67 21.69 10.32 13.63 6.22 method 1, case 2 V-0.962 V-0.960 V-0.951 V-0.936 V-0.907 V=0.852 V-0.737 N-900 N-450 N-150 N-60 N-30 N-40 N-40 CPU 425.34 214.73 74.74 32.25 30.56 87.87 169.94 Ax-0.10 method 2, case 2 V-0.965 V-0.966 V-0.968 not N-900 N-450 N-150 convergent CPU 584.5 429.57 284.78 method 3, case 2 V-0.964 V-0.964 V-0.961 V-0.955 V=0.945 over N-900 N-450 N-150 N-60 N-30 flow CPU 92.4 48.32 18.7 9.95 6.69 F method 1, case 3 V-0.962 V-0.959 V-O.952 V-0.936 V-0.906 V=0.853 V-0.777 N-900 N-450 N-150 N-60 N-30 N-40 N-40 CPU 678.98 341.2 117.27 49.53 46.19 138.48 269.18 Ax-0.063 method 2, case 3 V-0.965 V-0.965 V-0.965 not N-900 N-450 N-150 convergent CPU 792.4 498.34 286.5 method 3, case 3 V-0.963 V-0.963 V-0.960 V-0.954 V-0.944 over N-900 N-450 N-150 N-60 N—30 flow CPU 146.12 74.9 27.71 13.71 9.17 l. V-0.944 (analytical solution) 2. Method 1, analytical numerical method. 3. Method 2, Implicit method. 4. Method 3, Explicit the chemical reaction term and implicit for the convection diffusion term. Remarks: APPENDIX B PROGRAM LISTING FOR THE PILOT IGNITION MODEL ********************************************************************** COMPUTER MODEL FOR THE PILOT IGNITION OF WOOD by: Lin-Shyang Tzeng FIRDI in Taiwan P. 0. Box 246 Hsinchu 300,Taiwan R. 0. C. SEPTEMBER 1990 ********************************************************************** The following five lines is the example of input data (see next page under C "input data”) 0.0276 0.0000276 0.000276 0.00666 0.0037 0.00666 0.02 0 0 0 .218 1.0 1.657 0.866 9.59 18.2 3.3387 0.25 3.8 2.00 298. 0.6 0.197 6.087 .75 0.232 48000 560 110 21 221 241 3 40 10 0 241 20 1 0 0.135 1.657 0.69 0.07 0.01 31.0 PNM1F38.DAT The case shown above is in the file named PNM1F38.DAT combined analytical & numerical method solve chem reac. analytically when T(I). gt. 0.95 including the surface temp. solve the solid equation numerically not dumping the data after turn on the ignition source 132 133 Le number not - 1. PROGRAM ONEDPISN COMMON T(401),YF(401),YO(401),XX(401),A(6),B(6),D(6),N,NDUM l,NIGS,DTOD,IM,IN,IO,NODUF,TIME,DT,IFLAG,IDUMP,HL,NMAX,SMA l,I,DT1,DT2,DT3,CF,CO,AF,BT,Q,DA,FS,REAC,CLA,CSM,STP 1,AFU(6),DF(6),BF(6),YFS COMMON/WDATA/WTD(451),WTDO(451),WAMD(451),WAMDO(451),WDN(451) 1 ,WDNA(451),WDNC(451),WHCP(451),WCND(451),WA(451),WD(451),SIGM l ,WD1(451),WDNO(451),WB(451),WR(451),DNWF,DNF,HCPA,HCPC,CNDC 1 ,WDX,WTIME,EACT DIMENSION TOD(401),YFOD(401),YOD(401),TTD(401),YFTD(401) 1,YOTD(401) DOUBLE PRECISION WDN,WDNO,WDX CHARACTER*16 SFTPDAT OPEN(UNIT-10,FILE-'IGDUFLN.DAT',TYPE-‘NEW') OPEN(UNIT-1l,FILE-'DAISN.DAT',TYPE-'OLD') OPEN(UNIT-12,FILE-'OUPTN.DAT',TYPE-'NEW') OPEN(UNIT-l3,FILE-‘HTFXN.DAT',TYPE-'NEW') IGDUFLN.DAT data for check if there is ignition DAISN.DAT input data file 0UPTN.DAT file for temperature and species profile HTFXN.DAT file for wave of heat flux input data READ(11,*)DT,DT2,DT3,DX1,DX2,DX3,SM1,SM2,ST1,ST2,YFS,RLE READ(11,*)AF,BT,Q,DA,FS,HSF,CVF,TAMB,DNS,RKK,APS,EPSL,YOS READ(11,*)NMAX,NIG,IG,IM,IN,IO,IGN,NDUP,NHS,NW,IDUMP,NCH,M READ(11,*)DNF,HCPA,HCPC,CNDC,WDX,EACT 2400 2500 134 READ(11,2400)SFTPDAT FORMAT(A16) OPEN(UNIT-14,FILE-SFTPDAT,TYPE-'NEW') WRITE(14,*)DT,DT2,DT3,DX1,DX2,DX3,SM1,SM2,STl,ST2,YFS,RLE WRITE(14,*)AF,BT,Q,DA,FS,HSF,CVF,TAMB,DNS,RKK,APS,EPSL,YOS WRITE(14,*)NMAX,NIG,IG,IM,IN,IO,IGN,NDUP,NHS,NW,IDUMP,NCH,M WRITE(14,*)DNF,HCPA,HCPC,CNDC,WDX,EACT WRITE(14,2500)SFTPDAT FORMAT(lX,A16) DT - time step of no chem. react. DT2 - time step of premixed flame DT3 - time step of diffusion flame DXl - grid space near the surface DX2 - grid space in the middle area DX3 - grid space near the upper surface SM1 - const fuel flow rate SM2 - coeff. of linear fuel flow rate STl - constant surface temp. 3T2 - coeff. of linear surface temp. SM1,SM2,ST1,ST2 in this program are only dummy variable YFS - fuel concentration in the solid phase RLE - Lewis number AF - non-dimensional flame temperature BT - non-dimensional activation energy of fuel Q - non-dimensional heat of reaction of fuel DA - Damkohler number FS - stoichiometry ratio 135 HSF - radiation heat flux imposed on the wood CVF - convection heat transfer coeff. (W/m2 k) TAMB - ambient temperature k DNS - density of wood g/cm3 RKK - thermal conductivity of wood W/m k APS - pre-exponential factor of pyrolysis of wood EPSL - surface emissivity YOS - oxygen concentration of ambient air NMAX - maximum number of step NIG - number of time step to turn on the ignition source 16 - grid location of the ignition source IM - number of the grid near the surace IN - number of the grid in the middle area IO - total number of the grid IGN - number of grids for ignition source NDUP - number of step between calling subroutine DUMP NHS - number of step between turn on the ignition source NW - flag to control if write down data of temp. fuel and time IDUMP - number of data grids to be dumped NCH - interval of step to write down the surface temp. with respect to time. M - number of step to calculate the chemical reaction term DNF - final char density HCPA - specific heat of active wood HCPC - specific heat of char CNDC - thermal conductivity of char WDX - grid space of wood 136 EACT - activation energy of wood (kcal/g mole) DX1*(IM-1)+DX2*(IN-IM)+DX3*(IO-IN) - 1.00 for boundary layer thickness - 1.5 cm TRES-0.0011 TMAX-(DT+DT3)*0.5*NMAX+1.5 initial condition D0 100 I-l,IO T(I)-0. YF(I)-0. YO(I)-YOS 100 CONTINUE coefficient of matrix and boundary cond. B(1)-0. A(1)-(-1.)/DX1 AFU(1)--1./(DX1*RLE) NIGA-O NIGB-O NIGS-O NDUM-O NODUF-O NTD-O NTDl-O “APO NIGT-NIG DTTD-DT NDUPTD-NDUP NHSTD-NHS 102 104 106 110 108 137 IIGN-IG+(ICN/2) TIGPS-O DTOD-DT TIME-OT ERROl-0.0 ERRO-1.E-8 N-1 IFLAG-O INOWR-O specify the x-coordinate XX(I) (Variable grid) XX(1)-0. DO 102 I-2,IM XX(I)-XX(I-1)+DX1 CONTINUE DO 104 I-IM+1,IN XX(I)-XX(I-l)+DX2 CONTINUE DO 106 I-IN+1,IO XX(I)-XX(I-l)+DX3 CONTINUE send data to temporary file CONTINUE DO 108 I-l,IO TOD(I)-T(I) YFOD(I)-YF(I) YOD(I)-YO(I) CONTINUE C 111 138 start the new time step CONTINUE IF(TIME.GT.TMAX) GO TO 400 IF((MA-M).GT.2) GO TO 150 NCHl-(N/NCH)*NCH*NW NCH2-(N/NCH)*NCH IF(N.EQ.NCH1) WRITE(13,1100)TIME,T(2),T(IG) 1,T(IN),YF(2),YF(IG),YF(IN),Y0(2),YO(IG),YO(IN) 30 40 2200 CALL SFTN(SMA,STP,DT,TIME,HSF,CVF,TAMB,DNS,RKK,APS,EPSL,N) IF(N.EQ.NCH2) THEN locate the max temp IDUF-Z TDUF-O. DO 40 I-2,IO-1 IF(T(I).GT.TDUF) GO TO 30 GO TO 40 CONTINUE IDUF-I TDUF-T(I) CONTINUE TISEC-TIME/0.55187 WRITE(14,2200)TISEC,SMA,T(4),T(l),TDUF,IDUF,N END IF FORMAT(lX,'tsec,SMA,ST21,TPM,I,N-',F9.4,4F9.5,1X,I4,1X,15) IF(SMA.GT.2.) GO TO 400 T(l)-STP YF(l)-SMA*YFS 139 YO(l)-0. D(l)-SMA+(l./DX1) A(2)-DT*((SMA/DX1)-(1./(DX1**2))) D(2)-1.-(SMA*DT/DX1)+(2.*DT/(DX1**2)) B(2)-(-1.)*DT/(DX1**2) A(3)-SMA*DT/DX2-2.*DT/((DX1+DX2)*DX2) D(3)-1.+2.*DT/((DX1+DX2)*DX2)+2.*DT/((DX1+DX2)*DX1) 1-SMA*DT/DX2 B(3)-(-2.)*DT/((DX1+DX2)*DX1) A(4)-(SMA*DT/DX2)-DT/(DX2**2) D(4)-l.-(SMA*DT/DX2)+(2.*DT/(DX2**2)) B(4)-(-1.)*DT/(DX2**2) A(5)-SMA*DT/DX3-2.*DT/((DX2+DX3)*DX3) D(5)-1.+2.*DT/((DX2+DX3)*DX3)+2.*DT/((DX2+DX3)*DX2) 1-SMA*DT/DX3 B(5)-(-2.)*DT/((DX2+DX3)*DX2) A(6)-(SMA*DT/DX3)-DT/(DX3**2) D(6)-1.-(SMA*DT/DX3)+(2.*DT/(DX3**2)) B(6)-(-1.)*DT/(DX3**2) DF(1)-SMA+(1./(DX1*RLE)) AFU(2)-DT*((SMA/DX1)-(l./((DX1**2)*RLE))) DF(2)-1.-(SMA*DT/DX1)+(2.*DT/((DX1**2)*RLE)) BF(2)--1.*DT/((DX1**2)*RLE) AFU(3)-SMA*DT/DX2—2.*DT/(((DX1+DX2)*DX2)*RLE) DF(3)-l.+2.*DT/(((DX1+DX2)*DX2)*RLE)+2.*DT/(((DX1+DX2)*DX1)* 1 RLE)-SMA*DT/DX2 BF(3)--2 .*DT/( ( (DX1+DX2)*DX1)*RLE) 140 AFU(4)-(SMA*DT/DX2)-DT/((DX2**2)*RLE) DF(4)-1 . - (SMA*DT/DX2)+(2 . *DT/( (DX2**2)*RLE)) BF(4)--1.*DT/((DX2**2)*RLE) AFU(5)-SMA*DT/DX3-2.*DT/(((DX2+DX3)*DX3)*RLE) DF(5)-l . +2 . *DT/( ( (DX2+DX3)*DX3)*RLE)+2 . *DT/( ( (DX2+DX3) 1 *DX2)*RLE)-SMA*DT/DX3 BF(5)--2.*DT/(((DX2+DX3)*DX2)*RLE) AFU(6)-SMA*DT/DX3-DT/((DX3**2)*RLE) DF(6)-1.-(SMA*DT/DX3)+(2.*DT/((DX3**2)*RLE)) BF(6)--l.*DT/((DX3**2)*RLE) calculate the homogeneous chemical reaction DTl-DT/M MA-M D0 140 I-2,IO-1 assume that if T(I). less 1.E-12 ,then no chem reac. IF(T(I).LT.1.E-12) GO TO 140 IF(T(I).GT.0.95) GO TO 112 GO TO 113 112 CONTINUE CFl-Q*YF(I) C01-FS*Q*YO(I) CFO-CF1+COl CLA-O.5*(CFO+ABS(CF1-C01)) CSM-CFO-CLA CF-CF1+T(I) C0-C01+T(I) CALL AS 141 GO TO 121 113 CONTINUE CF-T(I)+(Q*YF(I)) CO-T(I)+(FS*Q*YO(I)) TB-T(I) IF(T(I).LT.0.2) GO TO 115 CALL CHEM T(I)-TB+REAC IF((CF-T(I)).LT.ERR0.0R.(CO-T(I)).LT.ERRO) THEN T(I)-TB GO TO 112 END IF GO TO 121 115 CONTINUE DO 120 J-1,MA TAFT(I) CALL CHEM AR-REAC T(I)-TA+(0.5*AR) CALL CHEM BRPREAC T(I)-TA+(0.5*BR) CALL CHEM CR-REAC T(I)-TA+CR CALL CHEM T(I)-TA+((AR+2.*BR+2.*CR+REAC)/6.) 142 IF((CF-T(I)).LT.ERROl.OR.(CO-T(I)).LT.ERROl) THEN T(I)-TB GO TO 112 END IF 120 CONTINUE 121 CONTINUE YF(I)-(CF-T(I) )/Q YO(I)-(C0-T(I) )/(FS*Q) 140 CONTINUE solve the temperature and species equation CALL SY IF(IFLAG.EQ.2) GO TO 150 GO TO 160 150 CONTINUE IF(INOWR.EQ.1) GO TO 152 WRITE(10,*)'IFLAG,N,DT—',IFLAG,N,DT 152 CONTINUE IFLAG-O DT-DT/2. NDUP-NDUP*2 NHS-NHS*2 NTDl-N+20 MA-O DO 155 I-l,IO T(I)-TOD(I) YF(I)-YFOD(I) YO(I)-YOD(I) 155 160 200 210 211 143 CONTINUE GO TO 110 CONTINUE check if there is ignition IF(N.LT.NIG) GO TO 200 GO TO 210 N—N+1 TIME-TIME+DT GO TO 110 CONTINUE IF(NIGS.EQ.1) GO TO 217 IF(N.LE.NIGT) GO TO 213 IF(T(IIGN).LT.0.15) GO TO 211 IF((TIME-TIMETD).GT.TRES.AND.T(IIGN).LT.0.2) GO TO 211 IF((N-NIGT).GT.4500.AND.T(IIGN).GT.0.4) GO TO 400 GO TO 213 not dump data, that pick the data at N-N+l CONTINUE NHS-NHSTD NIGT-N+NHS N-N+1 DT-DTTD NDUP-NDUPTD TIME-TIME+DT INOWR-O NCH-NCH/40 GO TO 110 144 213 CONTINUE IF((N-220).GT.NTD)DT-DT3 D0 214 I-2,IG IF(T(I).GT.1.0) THEN NIGS—l II-I WRITE(10,*)'pilot ign at n,t,i,T-’,N,TIME,I,T(I) TIMETDS-TIMETD-DTI WRITE(10,*)'NTD,TIMETD—',NTD,TIMETDS WRITE(10,*)'T,YF,YO,at t,I-',TTD(II),YFTD(II),YOTD(II) 1 ,TIMETDS,II WRITE(10,*)'T,YF,YO,at t,IG-',TTD(IG),YFTD(IC),YOTD(IC) l ,TIMETDS,IG WRITE(13,*)'pilot ign at n,t,i,T-',N,TIME,I,T(I) WRITE(13,*)'NTD,TIMETD-',NTD,TIMETDS WRITE(13,*)'T,YF,YO,at t,I-',TTD(II),YFTD(II),YOTD(II) 1 ,TIMETDS,II WRITE(13,*)'T,YF,YO,at t,IG-',TTD(IG),YFTD(IG),YOTD(IC) 1 ,TIMETDS,IG END IF 214 CONTINUE DO 216 I-IG+IGN+1,IN IF(T(I).GT.0.99) THEN NIGS-l II-I WRITE(10,*)'pilot ign at n,t,i,T-',N,TIME,I,T(I) TIMETDS-TIMETD-DTI 216 217 220 145 WRITE(10,*)’NTD,TIMETD-',NTD,TIMETDS WRITE(10,*)'T,YF,YO,at t,I=',TTD(II),YFTD(II),YOTD(II) ,TIMETDS,II WRITE(10,*)'T,YF,YO,at t,IG-',TTD(IG),YFTD(IG),YOTD(IC) ,TIMETDS,IG WRITE(13,*)'pilot ign at n,t,i,T-',N,TIME,I,T(I) WRITE(13,*)'NTD,TIMETD-',NTD,TIMETDS WRITE(13,*)'T,YF,YO,at t,I—',TTD(II),YFTD(II),YOTD(II) ,TIMETDS,II WRITE(13,*)'T,YF,YO,at t,IG-',TTD(IG),YFTD(IG),YOTD(IG) ,TIMETDS,IG END IF CONTINUE CONTINUE IF((N-400).GT.NTD) DT-DT3 NB-(N/NDUP)*NDUP IF(NB.EQ.N) CALL DUMP IF(NB.EQ.N) WRITE(10,*)'N,DT-',N,DT IF(NIGS.EQ.1) GO TO 232 turn on the ignition source IF(N.EQ.NIGT)GO TO 220 GO TO 232 CONTINUE DO 225 I-1,IO TTD(I)-T(I) YFTD(I)iYF(I) YOTD(I)-YO(I) 225 230 232 400 146 CONTINUE NTD-N+1 NCH-NCH*40 TIMETD-DT+TIME if TDUF gt 0.28 not turn on ign source DT-DT2 DO 230 I-1,IGN .T(IG+I)-1. CONTINUE CONTINUE N-N+1 TIME-TIME+DT IF(N.GT.NMAX) GO TO 400 IF(N.EQ.NTD1) THEN DT-DT2 NDUP-NDUPTD NHS-NHSTD INOWRPI END IF GO TO 110 CONTINUE REWIND(UNIT-12) WRITE(10,*)'N,t,DT-',N,TIME,DT WRITE(12,*)'data for t,N-',TIMETD,NTD D0 410 I-2,IO WRITE(12,1200)I,XX(I),TTD(I),YFTD(I),YOTD(I) 410 CONTINUE 147 1100 FORMAT(' t,TF02,IG,IN',F8.3,1X,9F5.2) 1200 FORMAT(' I,X,T,YF,YO=',IS,lX,Fll.4,3F7.4) NN-NMAX+450 WRITE(12,*)' NMAX+450-',NN WRITE(13,*)' NMAX+450-',NN IF(IFLAG.EQ.1) WRITE(10,*)'IFLAG,N-',IFLAG,N CLOSE(UNIT-lO) CLOSE(UNIT-ll) CLOSE(UNIT-12) CLOSE(UNIT-13) STOP END ************************************** This subroutine is to calculate the solution of tri-diagonal matrix SUBROUTINE SY ************************************** COMMON T(401),YF(401),YO(401),XX(401),A(6),B(6),D(6),N,NDUM l ,NIGS,DTOD,IM,IN,IO,NODUF,TIME,DT,IFLAG,IDUMP,HL,NMAX,SMA l ,I,DT1,DT2,DT3,CF,CO,AF,BT,Q,DA,FS,REAC,CLA,CSM,STP 1 ,AFU(6),DF(6),BF(6),YFS DIMENSION DD(401),DD1(401),AA(401),BB(401) D0 10 I-2,IM-l DD(I)-D(2) AA(I)-A(2) BB(I)-B(2) 10 CONTINUE DD(IM)-D(3) 20 30 148 AA(IM)-A(3) BB(IM)-B(3) DO 20 I-IM+1,IN-1 DD(I)-D(4) AA(I)-A(4) BB(I)-B(4) CONTINUE DD(IN)-D(5) AA(IN)-A(5) BB(IN)-B(5) DO 30 I-IN+1,IO-l DD(I)-D(6) AA(I)-A(6) BB(I)-B(6) CONTINUE solve the temp equ. DD(1)-1. AA(1)-0. DD(IO)-l. BB(IO)-0. solve by Gauss elimination DD1(1)-DD(1) DO 40 I-2,IO-l CC-BB(I)/DD1(I-1) DD1(I)-DD(I)-(CC*AA(I-1)) T(I)-T(I)-(CC*T(I-l)) 40 CONTINUE 50 52 54 149 DO 50 I-2,IO-1 II-IO-I+1 T(II)-(T(II)-AA(II)*T(II+1))/DD1(II) IF(T(II).LT.1.E-22)T(II)=O. CONTINUE solve the fuel & Oxygen equ. DD(1)-DF(1) AA(1)—AFU(1) DD1(1)-DD(1) DO 52 I-2,IM-1 DD(I)-DF(2) AA(I)-AFU(2) BB(I)-BF(2) CONTINUE DD(IM)-DF(3) AA(IM)-AFU(3) BB(IM)-BF(3) DO 54 I-IM+1,IN-1 DD(I)-DF(4) AA(I)-AFU(4) BB(I)-BF(4) CONTINUE DD(IN)-DF(5) AA(IN)-AFU(5) BB(IN)-BF(5) DO 56 I-IN+1,IO-1 DD(I)-DF(6) 150 AA(I)-AFU(6) BB(I)-BF(6) S6 CONTINUE DO 60 I-2,IO-1 CC-BB(I)/DD1(I-l) DD1(I)-DD(I)-(CC*AA(I-l)) YF(I)-YF(I)-(CC*YF(I-1)) YO(I)-YO(I)-(CC*YO(I-l)) 60 CONTINUE DO 70 I-2,IO II-IO-I+l YF(II)-(YF(II)-AA(II)*YF(II+1))/DD1(II) YO(II)-(YO(II)-AA(II)*YO(II+1))/DD1(II) IF(YF(II).LT.1.E-22)YF(II)-0. IF(YO(II).LT.1.E-22)YO(II)-0. 70 CONTINUE RETURN END ****************************************** This subroutine is to dump the data SUBROUTINE DUMP ****************************************** COMMON T(401),YF(401),Y0(401),XX(401),A(6),B(6),D(6),N,NDUM 1 ,NIGS,DTOD,IM,IN,IO,NODUF,TIME,DT,IFLAG,IDUMP,HL,NMAX,SMA 1 ,I,DT1,DT2,DT3,CF,CO,AF,BT,Q,DA,FS,REAC,CLA,CSM,STP l ,AFU(6),DF(6),BF(6),YFS IA-IO/IDUMP 151 HX-0.8299 NDUMFNDUM+1 IF(NDUM.GT.1.5) GO TO 20 IF(NIGS.EQ.1) GO TO 5 NDUM-O REWIND(UNIT=12) 5 CONTINUE IF(NIGS.EQ.1) WRITE(10,*)'Pilot ignition at N,T-',N,TIME IF(NIGS.EQ.1) WRITE(13,*)'Pilot ignition at N,T—',N,TIME WRITE(12,*) ' N,TIME,DT-',N,TIME,DT DO 10 I-2,IO,IA HFX-HX*(T(I)-T(I-1))/(XX(I)-XX(I-1)) WRITE(12,2000)I,XX(I),T(I),YF(I),YO(I),HFX 10 CONTINUE GO TO 70 20 CONTINUE IF(N.EQ.NMAX) GO TO 5 check if there is diffusion flame IDUF-2 TDUF-O. DO 40 I-2,IO-1 IF(T(I).GT.TDUF) GO TO 30 GO TO 40 30 CONTINUE IDUF-I TDUF-T(I) 4O CONTINUE 60 152 IF(TDUF.LT.0.28) THEN WRITE(10,*)'f1ame quenched at N,T-',N,TIME WRITE(13,*)'flame quenched at N,T-',N,TIME NDUMFO NIGS-O NODUF-O NMAX-N+25000 REWIND(UNIT=12) GO TO 5 END IF IF(((NDUM/20)*20).EQ.NDUM) GO TO 60 IF(NODUF.EQ.1) GO TO 70 NODUF-O CC-O. IF(TDUF.GT.0.6.AND.YO(1).LT.0.001) THEN DT-DT3 NMAX~N+4500 NODUF-l END IF CONTINUE IF(NODUF.EQ.0) GO TO 70 WRITE(10,*)'There is a diffusion flame at' WRITE(10,*)'I,N,TIME-',IDUF,N,TIME WRITE(10,*)'SMA,STP-',SMA,STP WRITE(13,*)'There is a diffusion flame at' WRITE(13,*) 'I,N,TIME,DT-',IDUF,N,TIME,DT,'HX, J/CM2 S' I-IDUF 65 70 2000 153 HFX-HX*(T(I) -T(I-1))/(XX(I) -XX(I-1)) WRITE(13,2000)I,XX(I),T(I),YF(I),YO(I),HFX WRITE(12,*) ' N,TIME,DT-',N,TIME,DT DO 65 I-2,IO,IA HFX-HX*(T(I)-T(I-1))/(XX(I)-XX(I-1)) WRITE(12,2000)I,XX(I),T(I),YF(I),YO(I),HFX CONTINUE CONTINUE FORMAT(' I,X,T,YF,YO-',IS,1X,F11.4,3F7.4,F9.4) RETURN END **************************************** This subroutine is to calculate the chemical reaction rate SUBROUTINE CHEM mmm*******m*w*mm****** COMMON T(401),YF(401),YO(401),XX(401),A(6),B(6),D(6),N,NDUM l ,NIGS,DTOD,IM,IN,IO,NODUF,TIME,DT,IFLAG,IDUMP,HL,NMAX,SMA 1 ,I,DTl,DT2,DT3,CF,CO,AF,BT,Q,DA,FS,REAC,CLA,CSM,STP l ,AFU(6),DF(6),BF(6),YFS C1-1.-T(I) IF(T(I).LT.1.E-2O) THEN CO TO 112 END IF C2-BT*Cl/(l.-AF*C1) IF(ABS(CZ).GT.30.) CO TO 112 C3-1./EXP(ABS(CZ)) IF(C2.LT.O.) C3-l./C3 154 CO TO 114 112 C3-0. IF(CZ.LT.0.) IFLAG-l 114 CONTINUE C4-C3*DT1/(FS*Q) REAC-DA*(CF-T(I))*(CO-T(I))*C4 RETURN END ******************************** This subroutine is to calculate the chemical reaction rate near flame temperature SUBROUTINE AS ******************************** COMMON T(401),YF(401),YO(401),XX(401),A(6),B(6),D(6),N,NDUM l ,NIGS,DTOD,IM,IN,IO,NODUF,TIME,DT,IFLAG,IDUMP,HL,NMAX,SMA 1 ,I,DT1,DT2,DT3,CF,CO,AF,BT,Q,DA,FS,REAC,CLA,CSM,STP 1 ,AFU(6),DF(6),BF(6),YFS Cl-l.-T(I) IF(T(I).LT.1.E-20) THEN GO TO 112 END IF C2-BT*Cl/(l.-AF*Cl) IF(ABS(CZ).GT.30.) GO TO 112 C3-1./EXP(ABS(C2)) IF(C2.LT.0.) C3-1./C3 GO TO 114 112 C3-0. 155 IF(C2.LT.0.) IFLAG=1 114 CONTINUE C4-DA*CLA*DT/(FS*Q) IF(C4.LT.1.) GO TO 116 IF(C3.GT.(50./C4)) THEN T(I)-T(I)+CSM GO TO 120 END IF 116 CONTINUE C5-C4*C3 C6-1./EXP(C5) T(I)-T(I)+CSM*(l.-C6) 120 CONTINUE RETURN END C this subroutine solves the solid equations numerically SUBROUTINE SFTN(SMA,STP,DT,TIME,HSF,CVF,TAMB,DNW,CNDA,APS,EPSL,N) c ***************************************************************** COMMON/WDATA/WTD(451),WTDO(451),WAMD(451),WAMDO(451),WDN(451) 1 ,WDNA(451),WDNC(451),WHCP(451),WCND(451),WA(451),WD(451),SIGM l ,WD1(451),WDNO(451),WB(451),WR(451),DNWF,DNF,HCPA,HCPC,CNDC 1 ,WDX,WTIME,EACT DOUBLE PRECISION CCl,CC2,CC3,WDT,WDX,WDN,WDNO C DNW- Virgin wood density C DNFb Final char density C HCPA- Specific heat active wood C HCPC- Specific heat of char 156 CNDA- Thermal condcutivity of active wood CNDC- Thermal conductivity of char EPSL- Emissivity HSF- Imposed heat flux CVFB h/EPSL ,h is convective heat transfer coefficient TAMB- Ambient temperature SIGN? 5.6696E-12 APS- Ad, pre-exponential factor of wood NMAXP Maximum.number of step WDXP Grid space of wood DT- Time step IO- number of grid space initial condition IO-350 IF(N.GT.1) GO TO 200 change the unit of input data CNDAF0.01*CNDA CNDC-0.01*CNDC CVF-0.0001*CVF/EPSL EACT-EACT*1000./1.987 DO 100 I-1,IO WTD(I)-TAMB WTDO(I)-TAMB WDN(I)-DNW WDNO(I)-DNW WAMD(I)-0. WAMDO(I)-0. 100 160 170 200 157 CONTINUE SICM-5.6696E—12 wnxz-wbx**2 DNWF-DNW-DNF CONTINUE WTIME-TIME DO 170 I-1,IO WDN(I)-WDNO(I) WTD(I)-WTDO(I) CONTINUE CONTINUE WDT-DT/0.55187 IF((WTIME-TIME).GT.1.E-12) CO TO 160 WDNA(1)-DNW*(WDN(1)-DNF)/DNWF WDNC(1)-DNF*(DNW-WDN(1))/DNWF WHCP(l)-(WDNA(1)*HCPA+WDNC(1)*HCPC)/WDN(1) WCND(1)-(WDNA(l)*CNDA/DNW)+(WDNC(1)*CNDC/DNF) WA(1)--WCND(l)/WDX2 WB(1)-WA(1) WD(1)-(-WA(l)-WB(1))+(WDN(1)*WHCP(1)/WDT) WR(1)-WDN(1)*WHCP(1)*WTD(1)/WDT DO 210 I-2,IO WDNA(I)-DNW*(WDN(I)-DNF)/DNWF WDNC(I)-DNF*(DNW-WDN(I))/DNWF WHCP(I)-(WDNA(I)*HCPA+WDNC(I)*HCPC)/WDN(I) WCND(I)-(WDNA(I)*CNDA/DNW)+(WDNC(I)*CNDC/DNF) WA(I)--WCND(I)/WDX2 158 WB(I)-WA(I-l) WD(I)-(-WA(I)-WB(1))+(WDN(I)*WHCP(I)/WDT) WR(I)-WDN(I)*WHCP(I)*WTD(I)/WDT 210 CONTINUE WA(1)-l. WD(1)--l. WB(1)-0. WA(IO)-0. WD(IO)-l. WB(IO)--1. WR(1)-(-WDX*EPSL/WCND(1))*(HSF-CVF*(WTD(1)-TAMB)-SIGM*((WTD(1)**4) 1 -(TAMB**4))) where CVF-h/EPSL WR(IO)-0. solve the temperature equation DO 220 I-1,IO WDl(I)-WD(I) 220 CONTINUE DO 230 I-2,IO CC—WB(I)/WD1(I-1) WDl(I)-WD1(I)-(CC*WA(I-1)) WR(I)-WR(I)-(CC*WR(I-l)) 230 CONTINUE WR(IO)-WR(IO)/WD1(IO) WTDO(IO)-WTD(IO) WTD(IO)-WR(IO) DO 240 I-2,IO 159 II-IO-I+l WR(II)-(WR(II)-WA(II)*WR(II+1))/WD1(II) WTDO(II)-WTD(II) WTD(II)-WR(II) 240 CONTINUE WTD(IO+1)-WTD(IO) D0 250 I-1,IO CC-EACT/WTD(I) CCl-APS/EXP(CC) APS-Ad CC2-CC1+(l./WDT) CC3-(WDN(I)/WDT)+CCl*DNF WDNO(I)-WDN(I) WDN(I)-CC3/CC2 250 CONTINUE DO 260 I-2,IO II—IO-I+1 WAMD(II)iWAMD(II+1)-WDX*(WDN(II)-WDNO(II))/WDT 260 CONTINUE STP-(WTD(1)-TAMB)/(2230.-TAMB) SMA-WAMD(1)*3230. WTIME-TIME RETURN END ***********‘k*************************‘k*********************~k*********** APPENDIX C PROGRAM LISTING FOR THE STEADY STATE DIFFUSION FLAME (FINITE DIFFERENCE METHOD AND ACTIVATION ENERGY ASYMPTOTICS) C.l Finite difference numerical code C PROGRAM FOR THE STEADY STATE DIFFUSION FLAME C by: Lin-Shyang Tzeng C solve T,YF,YO separately C for steady state diffusing flame C example for the input data of this program C 0.0037 0.0037 0.0037 0.2 0.1656 0.7 0 0.002 0 0 1.0 1.0 C 0.866 9.6 18.2 3.33E7 0.25 0.232 C 2000 2 101 181 191 273 0 1.0 10 l 10 10 C The result is writing on the file 'IGDUFLS.DAT' C In the ICDUFLS.DAT, the last line SMA,SMAO are the fuel flow rate , C SMA is the fuel flow rate where there is not a diffusion flame C SMAO is the fuel flow rate where there is a diffusion flame PROGRAM ONEDSL COMMON T(401),YF(401),YO(401),XX(401),A(6),B(6),D(6),N,SMA l ,RR(401),BT,AF,DA,FS,IM,IN,IO,NODUF,DX1,DX2,DX3,IFLAG,IDUMP 2,PM1,NDST,YFS,RLE,AFU(6),BF(6),DF(6),ERR COMMON DDF(401),AAF(401),BBF(401),DD(401),DD1(401),AA(401) l ,BB(401),DDO(401),AAO(401),BBO(401) 160 161 DIMENSION TOD(401),YFOD(401),YOD(401) OPEN(UNIT-10,FILE='IGDUFLS.DAT',TYPE-'NEW') OPEN(UNIT-11,FILE-‘DAINS.DAT',TYPE—'OLD') OPEN(UNIT-12,FILE='OUPTS.DAT',TYPE-‘NEW') input data READ(11,*)DX1,DX2,DX3,SMA,STP,PARM,PM1,PSM,IPM,MX,YFS,RLE READ(11,*)AF,BT,Q,DA,FS,YOS READ(11,*)NMAX,NDST,IG,IM,IN,IO,IRA,RAD,NHS,NW,IDUMP,NCH IRA-l including the radiation NDST-l the density and the conductivity is not constant RLE ; Levise number FARM ; the ratio of the old and new data for the iteration PMl ; maximum error for convergerne PSM ; error for the fuel flow rate of two NMAX ; maximum number for iteration IG ; initial flame location IDUMP ; number of dumped data (IA-IO/IDUMP) RAD Planck mean absorption coeff Kp (/m) NCH,NHS undefined HL-(IM-1)*DX1+(IN-IM)*DX2+(IO-IN)*DX3 HL2-HL**2 XF-l.-(LOC(0.058/YFS+1.))/SMA IG-IO*XF initial condition SMAO-SMA SMB-O. ERO-l.E-3 102 104 106 107 162 coefficient of matrix and boundary cond. B(1)-0. A(1)-(-1.)/DX1 AFU(1)-A(1)/RLE N-l NL-O IFLAG-O XX(1)-0. DO 102 I-2,IM XX(I)-XX(I-1)+DX1 CONTINUE DO 104 I-IM+1,IN XX(I)-XX(I-1)+DX2 CONTINUE DO 106 I-IN+1,IO XX(I)-XX(I-1)+DX3 CONTINUE initial condition DO 107 I-1,IG-3 T(I)-XX(I)/XX(IG) YF(I)-SMA*(1.-T(I)) YO(I)-0. T(I)-T(I)+STP CONTINUE DO 108 I-IG+3,IO m>--mw«ma» YF(I)-0. 108 109 110 111 163 YO(I)-YOS*(1.-T(I)) CONTINUE DO 109 I-l,7 II-I-3 T(IC+II)-l. YF(IG+II)-YF(IG-4) YO(IG+II)-YO(IG+4) CONTINUE start the new iteration CONTINUE IF(SMA.LT.1.E-9) GO TO 400 D0 111 I-1,IO TOD(I)-T(I) YFOD(I)-YF(I) YOD(I)-Y0(I) CONTINUE IF(IPM.EQ.1) GO TO 400 ERR-0. calculate the old chemical reaction CALL CHEM DO 120 I-2,IO-1 RTP-O. IF(IRA.EQ.1)THEN TP-T(I)*l932.+298 DNST-0.3528/TP RTP-RAD*(TP**4)*1.576E-l6/DNST END IF 120 130 135 164 RR(I)-Q*YF(I)*YO(I)*RR(I)-RTP IF(RR(I).LT.0.) RR(I)-O. CONTINUE CALL COEF IF(IFLAG.EQ.1) GO TO 400 solve the diffusion CALL SY D0 130 I-2,IO-l ERFABS(T(I)'RR(I)) IF(ER.GT.ERR) ERR-ER T(I)-(1.-PARM)*T(I)+PARM*RR(I) CONTINUE CALL CHEM CALL COEF IF(((N/2)*2).EQ.N) GO TO 142 CALL SYF DO 135 I-l,IO-1 ER-ABS(YF(I)-RR(I)) IF(ER.GT.ERR) ERR-ER YF(I)éYF(I)*(l.-PARM)+PARM*RR(I) CONTINUE CALL CHEM CALL COEF CALL SYO D0 140 I-1,IO-1 ER-ABS(Y0(I)-RR(I)) IF(ER.GT.ERR) ERR-ER 165 YO(I)-YO(I)*(1.-PARM)+PARM*RR(I) 140 CONTINUE GO TO 148 142 CONTINUE CALL SYO DO 145 I-1,IO-1 ER-ABS(YO(I)-RR(I)) IF(ER.GT.ERR) ERR-ER YO(I)-YO(I)*(1.-PARM)+PARM*RR(I) 145 CONTINUE CALL CHEM CALL COEF CALL SYF D0 147 I-l,IO-l ER-ABS(YF(I)-RR(I)) IF(ER.GT.ERR) ERR-ER YF(I)-YF(I)*(1.-PARM)+PARM*RR(I) 147 CONTINUE 148 CONTINUE IF(IFLAG.EQ.2) ERR-l. N-N+1 IF(N.GT.NMAX) THEN WRITE(10,*)'not converge for N=’,N GO TO 400 END IF IF(ERR.GT.ERO) GO TO 111 BRO-1.E-4 210 220 166 N-l TMX-O. II-o DO 210 I-2,IO-l IF(T(I).GT.TMX) THEN TMX-T(I) II-I END IF CONTINUE WRITE(10,*)'SMA,TMX,II,X-',SMA,TMX,II,XX(II) IF(NW.EQ.0) CO TO 410 IF(MX-l)220,230,240 CONTINUE IF(TMX.GT.0.27) THEN SMAO-SMA SMA-0.5*(SMA+SMB) SMR-ABS(l.-(SMA/SMAO)) IF(SMR.LT.PSM) THEN SMA-SMA*0.5 SMB-SMA END IF N-l GO TO 110 END IF DO 225 I-l,IO T(I)-TOD(I) YF(I)dYFOD(I) 225 230 232 234 236 167 YO(I)-YOD(I) CONTINUE SMB-SMA SMA-0.5*(SMAO+SMA) SMR-ABS(l.-(SMA/SMAO)) IF(SMR.LT.PSM) GO TO 400 N-l GO TO 111 CONTINUE IF(TMX.GT.0.1) THEN DXl-DX1-0.003S Dx2-Dx2-o.00035 Dx3-Dx3-o.oo35 A(1)-(21.)/DX1 DO 232 I-2,IM XX(I)-XX(I-1)+DX1 CONTINUE DO 234 I-IM+1,IN XX(I)-XX(I-1)+DX2 CONTINUE DO 236 I-IN+1,IO XX(I)-XX(I-1)+DX3 CONTINUE N-I GO TO 110 END IF GO TO 400 168 240 CONTINUE WRITE(10,*)'T(1),YF(1)=',T(1),YF(1) IF(TMX.GT.0.27) THEN YFS-YFS-0.01 N-l GO TO 110 END IF 400 CONTINUE SMA-SMA*1000.0/3230.0 SMAO-SMAO*1000.0/3230.0 «IF(MX.EQ.0) WRITE(10,*)'minimum fuel flow rate (mg/cm2 sec),SMA l ,SMAO-',SMA,SMAO IF(MX.EQ.1) WRITE(10,*)'minimum quench dist.XX(IO),DX1,DX2 1 -',XX(IO),DX1,DX2 . IF(MX.EQ.2) WRITE(10,*)'min fuel con YFS=',YFS IF(IFLAG.EQ.2) GO TO 410 IF(ERR.GT.0.01) THEN WRITE(10,*)'ERR-',ERR GO TO 410 END IF DO 410 I-2,IO-1 T(I)-TOD(I) YF(I)-YFOD(I) Y0(I)-YOD(I) 410 CONTINUE CALL DUMP IF(IFLAG.EQ.2) WRITE(10,*)'IFLAG,N-',IFLAG,N 169 IF(IFLAG.EQ.1) WRITE(10,*)'IFLAG,N=',IFLAG,N CLOSE(UNIT-lO) CLOSE(UNIT-ll) CLOSE(UNIT-l2) STOP END ************************************** This subroutine solve the tri-diagonal matrix SUBROUTINE SY ************************************** COMMON T(401),YF(401),YO(401),XX(401),A(6),B(6),D(6),N,SMA l ,RR(401),BT,AF,DA,FS,IM,IN,IO,NODUF,DX1,DX2,DX3,IFLAG,IDUMP 2 ,PM1,NDST,YFS,RLE,AFU(6),BF(6),DF(6),ERR COMMON DDF(401),AAF(401),BBF(401),DD(401),DD1(401),AA(401) l ,BB(401),DDO(401),AAO(401),BBO(401) RR(1)-T(1) RR(IO)-T(IO) IF(NDST.EQ.1) THEN DO 100 I-2,IM-l CDT-900./(T(I)*1932.+298.) CDA-900./(T(I+l)*1932.+298.) CXT-1./(DX1*DX1*CDT) CXAP1./(DX1*DX1*CDA) SMX-SMA/DXl DD(I)-CXA+CXT-SMX AA(I)-SMX-CXA BB(I)--CXT 100 200 170 CONTINUE CDT-900./(T(IM)*1932.+298.) CDA-900./(T(IM+1)*1932.+298.) CXT-l./(DX2*DX1*CDT) CXA-l./(DX2*DX2*CDA) SMX-SMA/DXZ DD(IM)-CXA+CXT-SMX AA(IM)-SMX-CXA BB(IM)--CXT DO 200 I-IM+1,IN-l CDT~900./(T(I)*1932.+298.) CDA-900./(T(I+1)*1932.+298.) CXT-1./(DX2*DX2*CDT) CXA-1./(DX2*DX2*CDA) SMX-SMA/DXZ DD(I)-CXA+CXT-SMX AA(I)-SMX-CXA BB(I)--CXT CONTINUE CDT-900./(T(IN)*1932.+298.) CDA-900./(T(IN+1)*1932.+298.) CXT-1./(DX2*DX3*CDT) CXA-l./(DX3*DX3*CDA) SMX-SMA/DX3 DD(IN)-CXA+CXT-SMX AA(IN)-SMX-CXA BB(IN)--CXT 300 10 20 171 D0 300 I-IN+1,IO-1 CDT-900./(T(I)*l932.+298.) CDA-900./(T(I+1)*1932.+298.) CXT-1./(DX3*DX3*CDT) CXA-1./(DX3*DX3*CDA) SMXPSMA/DX3 DD(I)-CXA+CXT-SMX AA(I)-SMX-CXA BB(I)--CXT CONTINUE GO TO 35 END IF DO 10 I-2,IM-l DD(I)-D(2) AA(I)-A(2) BB(I)-B(2) CONTINUE DD(IM)-D(3) AA(IM)-A(3) BB(IM)-B(3) DO 20 I-IM+1,IN-1 DD(I)-D(4) AA(I)-A(4) BB(I)-B(4) CONTINUE DD(IN)-D(5) AA(IN)-A(5) 30 35 40 50 172 BB(IN)-B(5) DO 30 I-IN+1,IO-1 DD(I)-D(6) AA(I)-A(6) BB(I)-B(6) CONTINUE CONTINUE solve the temp equ. DD(1)-l. AA(l)-0. DD(IO)-1.. BB(IO)-0. solve by Gauss elimination DD1(1)-DD(1) DO 40 I-2,IO-l CC-BB(I)/DD1(I-1) DD1(I)-DD(I)-(CC*AA(I-1)) RR(I)-RR(I)-(CC*RR(I-1)) CONTINUE D0 50 I-2,IO-l II-IO-I+1 RR(II)-(RR(II)-AA(II)*RR(II+1))/DD1(II) IF(RR(II).LT.(-l.E-12)) THEN RR(II)-PM1*RR(II) IFLAG-Z END IF CONTINUE 80 10 30 173 GO TO 80 CONTINUE RETURN END ****************************************** This subroutine dump the data to a temporary file SUBROUTINE DUMP W‘ki’i’i’ifik******************************** COMMON T(401),YF(401),Y0(401),XX(401),A(6),B(6),D(6),N,SMA l ,RR(401),BT,AF,DA,FS,IM,IN,IO,NODUF,DX1,DX2,DX3,IFLAG,IDUMP 2 ,PM1,NDST,YFS,RLE,AFU(6),BF(6),DF(6),ERR COMMON DDF(401),AAF(401),BBF(401),DD(401),DD1(401),AA(401) 1 ,BB(401),DDO(401),AAO(401),BBO(401) IA-IO/IDUMP DO 10 I-l,IO,IA WRITE(12,2000)I,XX(I),T(I),YF(I),YO(I) CONTINUE check if there is diffusion flame NODUF-O IDUF-2 TDUF-O. DO 40 I-2,IO-l IF(T(I).GT.TDUF) GO TO 30 GO TO 40 CONTINUE IDUF-I TDUFhT(I) 174 40 CONTINUE WRITE(10,*)'max temp N,T,I-',N,TDUF,IDUF 2000 FORMAT(' I,X,T,YF,Y0=',15,1X,F11.4,3F7.4) RETURN END ************************************* This subroutine calculate the chemical reaction rate SUBROUTINE CHEM ************************************** COMMON T(401),YF(401),YO(401),XX(401),A(6),B(6),D(6),N,SMA 1 ,RR(401),BT,AF DA,FS,IM,IN,IO,NODUF,DX1,DX2,DX3,IFLAG,IDUMP 2 ,PM1,NDST,YFS,RLE,AFU(6),BF(6),DF(6),ERR COMMON DDF(401),AAF(401),BBF(401),DD(401),DD1(401),AA(401) 1 ,BB(401),DDO(401),AAO(401),BBO(401) D0 120 I-2,IO-l Cl-1.-T(I) C2-BT*Cl/(l.-AF*C1) IF(ABS(C2).GT.50.) GO TO 112 C3-1./EXP(ABS(C2)) IF(C2.LT.0.) C3-1./C3 GO TO 114 112 C3-0. IF(CZ.LT.0.) IFLAG-l IF(IFLAG.EQ.1) PRINT *,'I,N,T-',I,N,T(I) 114 CONTINUE IF(NDST.EQ.1) THEN Density is not constant 175 DAl-DA*943.4/(T(I)*l932.+298.) GO TO 115 END IF DAl-DA 115 CONTINUE RR(I)-C3*DA1 120 CONTINUE RETURN END **************************************** This subroutine calculate the coefficient of the tri-diagonal matrix SUBROUTINE COEF **************************************** COMMON T(401),YF(401),YO(401),XX(401),A(6),B(6),D(6),N,SMA 1 ,RR(401),BT,AF,DA,FS,IM,IN,IO,NODUF,DX1,DX2,DX3,IFLAG,IDUMP 2 ,PM1,NDST,YFS,RLE,AFU(6),BF(6),DF(6),ERR COMMON DDF(401),AAF(401),BBF(401),DD(401),DDl(401),AA(401) l ,BB(401),DDO(401),AAO(401),BBO(401) D(1)-SMA+(l./DX1) A(2)-1.0*((SMA/DX1)-(l./(DX1**2))) D(2)-0.-(SMA*1.0/DX1)+(2.*1.0/(DXl**2)) B(2)-(-1.)*1.0/(DX1**2) A(3)-SMA*1.0/DX2-2.*l.0/((DX1+DX2)*DX2) D(3)-0.+2.*1.0/((DX1+DX2)*DX2)+2.*1.0/((DX1+DX2)*DX1) l -SMA*1.0/DX2 B(3)-(-2.)*1.0/((DX1+DX2)*DX1) A(4)-(SMA*1.0/DX2)-l.0/(DX2**2) 176 D(4)-0.-(SMA*1.0/DX2)+(2.*l.0/(DX2**2)) B(4)-(-1.)*1.0/(DX2**2) A(5)-SMA*1.0/DX3-2.*l.0/((DX2+DX3)*DX3) D(5)-0.+2.*1.0/((DX2+DX3)*DX3)+2.*1.0/((DX2+DX3)*DX2) 1 -SMA*1.0/DX3 B(5)-(-2.)*l.0/((DX2+DX3)*DX2) A(6)-(SMA*1.0/DX3)-l.0/(DX3**2) D(6)-0.-(SMA*1.0/DX3)+(2.*1.0/(DX3**2)) B(6)-(-1.)*1.0/(DX3**2) DF(1)-SMA+(1./(DX1*RLE)) AFU(2)-(SMA/DX1)-(l./((DX1**2)*RLE)) DF(2)-(2./((DX1**2)*RLE))-(SMA/DXl) BF(2)-(-l.)/((DX1**2)*RLE) AFU(3)-SMA/DX2-2./((DX1+DX2)*DX2*RLE) DF(3)-2./((DX1+DX2)*DX2*RLE)+2./((DX1+DX2)*DX1 1 *RLE)-SMA/DX2 BF(3)-(-2.)/((DX1+DX2)*DX1*RLE) AFU(4)-(SMA/DX2)-1./((DX2**2)*RLE) DF(4)-2./((DX2**2)*RLE)-(SMA/DXZ) BF(4)-(-l.)/((DX2**2)*RLE) AFU(5)-SMA/DX3-2./((DX2+DX3)*DX3*RLE) DF(5)-2./((DX2+DX3)*DX3*RLE)+2./((DX2+DX3)*DX2 1 *RLE)-SMA/DX3 BF(5)-(-2.)/((DX2+DX3)*DX2*RLE) AFU(6)-SMA/DX3-l./((DX3**2)*RLE) DF(6)-2./((DX3**2)*RLE)-SMA/DX3 BF(6)-(-l.)/((DX3**2)*RLE) 100 177 RETURN END ********************************** This subroutine solves the tri-diagonal matrix for fuel SUBROUTINE SYF ********************************** COMMON T(401),YF(401),Y0(401),XX(401),A(6),B(6),D(6),N,SMA 1 ,RR(401),BT,AF,DA,FS,IM,IN,IO,NODUF,DX1,DX2,DX3,IFLAG,IDUMP 2 ,PM1,NDST,YFS,RLE,AFU(6),BF(6),DF(6),ERR COMMON DDF(401),AAF(401),BBF(401),DD(401),DD1(401),AA(401) l ,BB(401),DDO(401),AAO(401),BBO(401) RR(1)-SMA*YFS RR(IO)-YF(IO) IF(NDST.EQ.1) THEN D0 100 I-2,IM-l CDT-900./(T(I)*1932.+298.) CDA-900./(T(I+l)*1932.+298.) CXT-l./(DX1*DX1*CDT*RLE) CXA-1./(DX1*DX1*CDA*RLE) SMXPSMA/DXI DDF(1)-CXA+CXT-SMX AAF(I)-SMX—CXA BBF(I)--CXT CONTINUE CDT-900./(T(IM)*1932.+298.) CDAF900./(T(IM+1)*1932.+298.) CXT-l . /(DX2*DX1*CDT*RLE) 5 200 178 CXA-l./(DX2*DX2*CDA*RLE) SMX-SMA/DXZ DDF(IM)-CXA+CXT-SMX AAF(IM)—SMX-CXA BBF(IM)--CXT DO 200 I-IM+1,IN-1 CDT-900./(T(I)*1932.+298.) CDA-900./(T(I+l)*1932.+298.) CXT-1./(DX2*DX2*CDT*RLE) CXA-l./(DX2*DX2*CDA*RLE) SMX-SMA/DX2 DDF(I)-CXA+CXT-SMX AAF(I)-SMX-CXA BBF(I)--CXT CONTINUE CDT-900./(T(IN)*1932.+298.) CDA-900./(T(IN+1)*1932.+298.) CXT-1./(DX2*DX3*CDT*RLE) CXA-1./(DX3*DX3*CDA*RLE) SMX-SMA/DX3 DDF(IN)-CXA+CXT-SMX AAF(IN)-SMX-CXA BBF(IN)--CXT DO 300 I-IN+1,IO-1 CDT-900./(T(I)*l932.+298.) CDA-900./(T(I+1)*1932.+298.) CXT-l./(DX3*DX3*CDT*RLE) 300 310 10 179 CXA-1./(DX3*DX3*CDA*RLE) SMX-SMA/DXB DDF(I)-CXA+CXT-SMX AAF(I)-SMX-CXA BBF(I)--CXT CONTINUE DO 310 I-2,IO-1 DDF(I)-RR(I)*YO(I)+DDF(I) CONTINUE CDA-900./(T(2)*1932.+298.) DDF(1)-(1./(RLE*CDA*DX1))+SMA AAF(1)--l./(RLE*CDA*DX1) GO To 32 END IF DO 10 I-2,IM-1 DDF(I)-RR(I)*YO(I)+DF(2) AAF(I)-AFU(2) BBF(I)-BF(2) CONTINUE DDF(IM)-RR(I)*YO(I)+DF(3) AAF(IM)-AFU(3) BBF(IM)-BF(3) DO 20 I-IM+1,IN-l DDF(I)-RR(I)*YO(I)+DF(4) AAF(I)-AFU(4) BBF(I)-BF(4) 20 CONTINUE 30 32 35 180 DDF(IN)-RR(I)*YO(I)+DF(5) AAF(IN)-AFU(5) BBF(IN)-BF(5) DO 30 I-IN+1,IO-1 DDF(I)-RR(I)*YO(I)+DF(6) AAF(I)-AFU(6) BBF(I)-BF(6) CONTINUE solve the fuel equ. DDF(1)-DF(1) AAF(1)-AFU(1) CONTINUE D0 35 I-2,IO-1 RR(I)-0. CONTINUE DDF(IO)-1. BBF(IO)-0. solve by Gauss elimination DD1(1)-DDF(1) DO 40 I-2,IO—1 CC-BBF(I)/DD1(I-1) DD1(I)-DDF(I)-(CC*AAF(I-1)) RR(I)-RR(I)-(CC*RR(I-l)) 4O CONTINUE DO 50 I-2,IO II-IO-I+l RR(II)-(RR(II)-AAF(II)*RR(II+1))/DD1(II) IF(RR(II).LT.(-l.E-12)) THEN RR(II)-PM1*RR(II) IFLAG-2 END IF SO CONTINUE RETURN END 181 ******************************** This subroutine solves the tri-diagonal matrix for oxygen SUBROUTINE SYO ********************************t COMMON T(401),YF(401),YO(401),XX(401),A(6),B(6),D(6),N,SMA l ,RR(401),BT,AF,DA,FS,IM,IN,IO,NODUF,DX1,DX2,DX3,IFLAG,IDUMP 2 ,PM1,NDST,YFS,RLE,AFU(6),BF(6),DF(6),ERR COMMON DDF(401),AAF(401),BBF(401),DD(401),DD1(401),AA(401) l ,BB(401),DDO(401),AAO(401),BBO(401) RR(1)-0. RR(IO)-YO(IO) IF(NDST.EQ.1) THEN DO 100 I-2,IM-1 CDT-900./(T(I)*1932.+298.) CDA-900./(T(I+l)*1932.+298.) CXT—l./(DX1*DX1*CDT*RLE) CXAp1./(DX1*DX1*CDA*RLE) SMXPSMA/DXl DDO(I)-CXA+CXT-SMX AAO(I)-SMX-CXA 100 200 182 BBO(I)--CXT CONTINUE CDT-900./(T(IM)*1932.+298.) CDA-900./(T(IM+1)*1932.+298.) CXT-1./(DX2*DX1*CDT*RLE) CXA-l./(DX2*DX2*CDA*RLE) SMX-SMA/DXZ DDO(IM)-CXA+CXT-SMX AAO(IM)-SMX-CXA BBO(IM)--CXT DO 200 I-IM+1,IN-l CDT-900./(T(I)*1932.+298.) CDA-900./(T(I+l)*l932.+298.) CXT-1./(DX2*DX2*CDT*RLE) CXA-l./(DX2*DX2*CDA*RLE) SMX-SMA/DXZ DDO(I)-CXA+CXT-SMX AAO(I)-SMX-CXA BBO(I)--CXT CONTINUE CDT-900./(T(IN)*1932.+298.) CDA-900./(T(IN+1)*1932.+298.) CXT-1./(DX2*DX3*CDT*RLE) CXA-l./(DX3*DX3*CDA*RLE) SMX-SMA/DX3 DDO(IN)-CXA+CXT-SMX AAO(IN)-SMX-CXA 300 310 10 183 BBO(IN)--CXT DO 300 I-IN+1,IO-1 CDT-900./(T(I)*1932.+298.) CDA-900./(T(I+1)*1932.+298.) CXT-l./(DX3*DX3*CDT*RLE) CXA-1./(DX3*DX3*CDA*RLE) SMX-SMA/DX3 DDO(I)-CXA+CXT-SMX AAO(I)-SMX-CXA BBO(I)--CXT CONTINUE DO 310 I—2,IO-1 DDO(I)-DDO(I)+(RR(I)*YF(I)/FS) CONTINUE CDA-900./(T(2)*1932.+298.) DDO(1)-(1./(RLE*CDA*DX1))+SMA AAO(l)--l./(RLE*CDA*DX1) GO TO 32 END IF DO 10 I-2,IM-1 DDO(I)-(RR(I)*YF(I)/FS)+DF(2) AAO(I)-AFU(2) BBO(I)-BF(2) CONTINUE DDO(IM)-(RR(I)*YF(I)/FS)+DF(3) AAO(IM)-AFU(3) BBO(IM)-BF(3) 20 30 32 35 184 D0 20 I-IM+1,IN-l DDO(I)-(RR(I)*YF(I)/FS)+DF(4) AAO(I)-AFU(4) BBO(I)-BF(4) CONTINUE DDO(IN)-(RR(I)*YF(I)/FS)+DF(5) AAO(IN)-AFU(5) BBO(IN)-BF(S) DO 30 I-IN+1,10-1 DDO(I)-(RR(I)*YF(I)/FS)+DF(6) AAO(I)-AFU(6) BBO(I)-BF(6) CONTINUE solve the Oxygen equ. DDO(1)-DF(1) AAO(1)-AFU(1) CONTINUE DO 35 I-2,IO-1 RR(I)-0. CONTINUE DDO(IO)-1. BBO(IO)-0. solve by Gauss elimination DD1(1)-DDO(1) D0 40 I-2,IO-1 CC-BBO(I)/DD1(I-1) DD1(I)-DDO(I)-(CC*AAO(I-l)) C.2 185 RR(I)-RR(I)-(CC*RR(I-1)) 40 CONTINUE DO 50 I-2,IO II-IO-I+l RR(II)-(RR(II)-AAO(II)*RR(II+1))/DD1(II) IF(RR(II).LTu(-l.E-12)) THEN RR(II)-PM1*RR(II) IFLAG-Z END IF SO CONTINUE RETURN END ******************************************************************** Program for the calculation of extintion Damkohler number (Activation energy asymptotics method) This program is to calculate the minimum fuel flow rate by using the extinction Damkohler number derived by Linan example for the input data, assume the fuel is metnane SMA-0.0002, YOX-0.232, HCM-1.5, TL—640.0, YFS=0.25, ACTV-48.0, AEXF-2.424E15, PROGRAM EXTIN FS-0.25 100 186 FS is the stoichiometry ratio of the fuel N-O PRINT *,'input SMA(g/cm2 sec),HCM(cm),YOX,TL(k),YFS' READ(S,*)SMA,HCM,YOX,TL,YFS SMA is the fuel flow rate at the surface HCM is the boundary layer thickness YOX is the Oxygen mass fraction TL is the surface temperature, (absolute temperature k) YFS is the fuel mass fraction beneath the solid surface SMD-0.000001 PRINT *,'input activation energy, pre-exponedtial factor' PRINT *,'kcal/mole,cm3/g sec' READ(5,*)ACTV,AEXF ACTV is the activation energy of the gas fuel AEXF is the frequency factor of the gas fuel CONTINUE N-N+1 IF(N.GT.1000) THEN PRINT *,’ NOT CONVERGE AT N-',N PRINT *,'REDAE,REDA',REDAE,REDA GO TO 220 END IF IF(SMA.LT.1.E-12) GO TO 200 SMAD-SMA*3230. YFSO-YFS—(YFS+FS*YOX)/EXP(SMAD) TPND-35155.*YFSO STPL-TL/TPND 200 187 STPU-298./TPND DL-3230 . *SMA AREA-FS*YOX/YFSO BETA-STPL-STPU TACE-ACTV/(0.001983*TPND) TDEF-STPL+(ARFA-BETA)/(1.+ARFA) EBS-(TDEF**2)/TACE GAMA-l.-2.*(ARFA-BETA)/(l.+ARFA) ZETO-1.-(1./EXP(DL)) DAMKO-4.*0.000374*AEXF*YFSO*5.61E-8/(0.323*(SMA**2)) CN1-4.*DAMKO*(ZETO**2)*(EBS**3) CN2-((1.-(ZETO*ARFA/(l.+ARFA)))**2)*((1.+ARFA)**2) CN3-1./EXP(TACE/TDEF) REDA-CN1*CN3/CN2 CN4-1.-GAMA REDAE-2.7183*(CN4-(CN4**2)+0.26*(CN4**3)+0.055*(CN4**4)) IF(REDA.GT.REDAE) THEN SMA-SMA-SMD N—N+1 GO TO 100 END IF CONTINUE IF(((REDAE-REDA)/(REDAE+REDA)).GT.0.01) THEN SMA-SMA+SMD SMD-SMD*0.1 GO TO 100 END IF 188 220 CONTINUE PRINT * , ' SMA , HCM , YOX , REDA , REDAE- ' , SMA , HCM , YOX , REDA , REDAE TDEFD-TDEF*TPND IF(N.EQ.1) PRINT *,'SMA too sma11,re-test again' PRINT *,'diffusion flame temp.-',TDEFD STOP END ********************************************************************** Appendix D ESTIMATION OF SOLID PHASE PARAMETERS FROM PILOTED IGNITION EXPERIMENTAL DATA D-l- W A series of piloted ignition experiments have been conducted by Abu-Zaid (1988). A combustion wind tunnel was used to perform these experiments in a well-controlled atmosphere under external radiation. This tunnel consists of three main sections, the inlet section, the test section, and the exhaust section; it is schematically shown in Figure D-l. The crossectional area for air flow is 15cm x 8.5 cm, and the wood sample (14.5 cm long, 7.5 cm width and 3.75 cm high) is placed 30cm downstream of the leading edge. These wood samples were carefully instrumented by surface, bottom and in-depth thermocouples. During the experiments, the heaters were turned on at the desired radiant heat flux, and were allowed to warm up for about 15-20 minutes. The wood sample was then placed on the weighing table with its top surface flat along the tunnel base. A small natural gas pilot flame was -used downstream of the tunnel to initiate piloted ignition. The output of gas analyzers, thermocouples and the load cell was stored and processed by a microcomputer in real time. 189 190 m=C.r<.5.=Z<: uUZuAaamDP E .25.: _scuEtoexo 2: Lo .5520 2382.8 To museum 339‘ :9: Ee< So “ms..— U , :2: 4093.28 «CU Ou—TNOU .II! 202.5503 <53 .5 52.: . 6...: 4.1 «5:15 96 mgizs‘ «<0 .83 2:... e2: Shite: v 22: «z «o _o .5525: :2: ...—8:; 298228 396 x3. .52: / 3580:23F eeagezfl. _ :2: «music Gauge moses c2505» 25 :95 .. I 5253 coo; \ 95E.» “AT—.4.” /' k 53:: .5358 \ . agz 220m / / xz<._. E< 35mg“.— =O_= II. Axoamgzou 205: Ill. I I. 3p. a ._ I? #2 ”<0 silk agar—km: E