‘ ‘ ”4-”‘3i 1. $4“ ~ -, ' ' . mvéamnxw gut ' 1ft! (g- 7;- ‘- .' 4w {J‘ék‘ : 5g ‘ _ . ‘ .' ,lfira“ I * " .‘ ~ 4fi%u~§i§=fi ‘m'.r*‘--»‘"‘1‘-§E - ‘3. :- _La.u- -| u .u cc. “1“?" ‘.... dw'igéif «33%. T‘sifi-fwgiwgtu «- a - ., . .héxm‘; “33:1 “in? _ .Vg’img; "‘ “est-55: * "Jignyxxr-n 1, u"""1\z. ‘7* “‘13" 513-- .. . ._ . _ ‘ .. A , v... .. ‘ g . , A ‘%-7.a«'.—1r—x~»§ntsszé‘ag * ‘- u ‘ -_~ . . . ~ — ~ ..» M 'r . . rig: . , . _m~ -u hwy-.3. . f a; v ‘ . , ‘ ,. ‘ ‘5, ‘ 7' mt“:- MJ‘ .11.... “a; ‘ ‘ . ,. - .,. ¢ . ‘. may“. I' "H". HAW?“ - ‘4‘ , . A . 4 - ‘ 31.. «0 1.2.. . _ . "t:- 2, .23-,“ on.“ _ -' . , ‘ » . .55 .. -.~ «Y‘afiflwz-‘DA. - <"’vfl7"'.‘{‘¢‘3ltan 3 '3‘ > . . 3w v ~- fr‘.fii§$¢rufiz}‘tfa ' ' 4- 3:3? 3:41”: -., , 'p “N qwfi 4'55 ,4 1;; 4:7- , a 4‘ \ . 3-0 4‘ w __' . o . .H .~.~ .pev’ . . ‘4.» w — ‘ »' " . - ##me $7323: fa" «we 11¢ $9. ~. ; u. - fr}? 9' M§fi¥r«?? ”6%? 7 ,-':'._ ' v.93 .zg’L‘J-{VS;““‘. .:%W. “ ‘64? :. '5 air ‘7 ' m,“ h’ggudar $42k " w ' w ‘ ‘ ‘23‘ "3‘2“ 15" ' :1 u , If '2? 15’?“ g g; km” :2. u. 'A . _. «n y» ' l . -u -. . _ | ., 4‘ ‘ 715,1 1. .07 «fit 5'1 ‘ 5?- .ulfi": 4 v A'. ' - - u .‘ .A. - J.’ "31%;," ' . . v" flit J I . - _ _ ”"4 ans—fl- a :«é..."co3{_¢ ,3". 1'". ‘1‘ .'" .mlr‘ 3. ‘ ‘ u A .v- .. ’ ‘ ‘, :‘M141:‘ ,, “f" '5 ‘ h ‘ " . ~ ' ”:3 “'3‘ . "fit ~hut-n2 .' 'éifl" -‘ $5151.34 ,, ”3:4; ~.. ‘1 l. ‘.. 7‘ ‘- ‘fi 1: 3 in m . q: n];- {'1‘ H .t ,1“ Ufifiwfl‘n: . r3. '3 r .j~ [- ".' 1'“? .. 5"? 2' '31:") . ' l . ;.,$,,.!,tw In. \ ‘1 "- . ”71933 $031; '. 9 a: u, .flfiv'uni ‘ , may. :3 4‘" . :- n «a “ “‘67 ' : fi v1 5. I Z "5%?- flaw. T ”imam ll?lllllllml’illlfilL L 3 00917 5492 This is to certify that the dissertation entitled SELF-HEATING AND SELF-IGNITION IN DAIRY POWDERS presented by James Finbarr O ’ Connor has been accepted towards fulfillment of the requirements for Ph.D degree in Agricultural Engineering Date August 22, 1990 MS U is an Affirmative Action/Equal Opportunity Institution 0-12771 r LIBMKY Michigan State i University \ J —_ PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE H “W __ 211$ 1« 1 ifl" 3§% MSU Is An Affirmative Action/Equal Opportunity Imitation cAcirchmHJ SELF-HEATING AND SELF-IGNITION IN DAIRY POIDBRS by JAMES FINBARR O'CONNOR A DISSERTATION submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY in the Department of Agricultural Engineering 1990 ABSTRACT BY James Finbarr O’Connor Self-heating in dairy-based milk powders is a cause of fires and explosions in processing equipment and storage facilities in the dairy industry. Detection of fires in spray-dryers and siloes has not been successful. Prevention offers the best method of limiting the occurrence. Thus, it has become essential to predict the conditions which are conducive to the commencement and propagation of heating in milk powders. A numerical model was developed to simulate self-heating of a sphere of powder. The simulation takes into account variable environmental and product conditions which include finite surface and internal resistance to heat transfer. The model is solved by finite elements. The solution is accurate and stable. The finite element technique permits extension of the simulation to other geometries, including layers and irregularly-shaped particles of milk powder. The model was tested over the range of parameters and properties likely to be encountered commercially: values of Activation Energy from 40 kJ/kg mole to 100 kJ/kg mole and test temperatures from JamfitFnumerTCmmmm 127°C to 380'C were examined for sphere radii of 1.8 cm to 5.1 cm. The Differential Scanning Calorimeter (DSC) was employed to provide kinetic data for the simulation. The model was subsequently validated for predicting the Minimum Ignition Temperature (MIT) and the Time to Ignition (TtI) of the milk powder under the test conditions. In the experiments commercial powders were subjected to a high-temperature environment in a standard oven. Detailed time/temperature profiles of the spheres were obtained as they heated up to either an ignition or a non-ignition condition. Experimental inaccuracies in the DSC resulted in simulated MIT values averaging 13% above the experimental results. The model was further employed to estimate the MIT/TtI values using the oven-derived kinetic data. These simulated time/temperature profiles show excellent correlation with the experimental data. For instance, for a 2" sphere of skim-milk powder, the oven predicts an MIT of 161°C, compared to a predicted MODEL/OVEN value of 159'C, and a MODEL/DSC value of 179'C; hence the ’oven-based' model gives the best correlation with the experimental data. The model allows very precise prediction of the lowest -u 'r! I JamfianunerTCmmmm ignition temperature. It specifically identifies the lowest ambient temperature at which ignition occurs. This is a better combustion indicator than the MIT which is defined as the average of the ignition and non-ignition temperatures. Experimental and simulation results show that the TtI increases as the sphere radius increases. Also, the TtI decreases exponentially when the ambient temperature increases linearly. Below the MIT, ignition is not possible. The model predicts an exponential decrease in the TtI for a linear decrease in the Activation Energy. Above an Activation Energy value of 90 kJ/mol, no ignition is possible under the standard test conditions. If the surface heat transfer coefficient is 15 W/m2 K or greater, -ignition takes place at approximately 150 mins for a 4" sphere. As a result of the findings, a number of specific practical recommendations are made regarding the prevention of self-heating and self-ignition in dairy E. W. Bakker-Ark Major Professori%zayé27 ”WM Approved R.D. von Bernuth f/za/yo Chairperson powders. To : Marion, Aoife and Orla ACKNO'IBDGBHBNTS The author wishes to express his thanks to the many people who assisted him in the course of this work. In particular he would like to thank University College Cork and the Kellogg Corporation of Battle Creek, Michigan for affording him the opportunity to undertake this course of study at Michigan State University. He is particularly grateful to the members of his advisory committee for their sustained interest and support over the duration of the project: Professor Fred Bakker-Arkema who provided invaluable guidance, friendship and support throughout: Professor Ian Gray similarly provided inestimable support from the moment of our first arrival in Lansing. Dr. Dennis Heldman, who provided the initial impetus for the program and Professors Eric Grulke and Jim Steffe who closely followed the progress of the work and provided support as required. On the research end he is indebted to Professor Chris Synnott and Mr. Diarmaid MacCarthy of UCC's Food Engineering Department and Dr. Joe Buckley of the Food Technology Department and to Dr. Thomas Duane, Mr. Liam O'Donnell and. Ms. Aine Curtain for’ many useful discussions; Messrs. John Barrett, Denis Ring and Jeseph O'Mahony who provided valued technical assistance and Ms. Rita Kelleher for her secretarial help. The support of the UCC Computer Bureau staff is also acknowledged. Finally he is extremely grateful to his wife, Marion, for her unstinting support and assistance over the long hours of sometimes seemingly endless toil and never-ending miles that have brought this project to fruition. List of Tables x List of Figures xii 1. Introduction 1 1.1 The Cost to Industry 10 1.2 Thermodynamics 12 1.3 'Classical' Ignition Theory 15 1.3.1 Minimum Ignition Temperature 17 1.3.2 The Frank-Kamenetskii Dimensionless 17 Reaction Rate 2. Objectives 20 3. Literature Review 22 3.1 Traditional Fire/Combustion Theory 22 3.1.1 The Semenov Model 22 3.1.2 The Frank-Kamenetskii Model 26 3.1.2.1 Summary of Frank-Kamenetskii 31 Assumptions 3.1.2.2 Use of Frank-Kamenetskii Theory 31 3.1.3 The Validity of the Traditional Approximations 36 3.1.3.1 The Frank-Kamenetskii Critical Parameter 37 3.1.3.2 Modifications to the Basic F-K Model 40 3.2 Shape and other Factors Influencing Criticality 43 3.2.1 Varying Thermal Conductivity 46 3.2.2 Reactant Consumption 48 3.3 Self-ignition in Products other than Milk Powders 50 3.3.1 Self-ignition in Wool 52 3.4 Milk Powder Ignition Data 57 3.4.1 IDF Summary (1987) 62 3.5 Calorimetry and Kinetic Studies 62 3.5.1 Reaction Rates 63 3.5.2 Calorimetry 67 3.5.3 DSC Calculation Procedure 69 3.5.4 DSC Analysis 73 3.5.5 Specific Heat 76 3.5.6 Density 77 3.6 Basic Mathematical Theory 78 shah-lb. eeee fihfibhhbbhfibbhbafihh NNH MUUIMUIUIUIUIUIUIUIMUIMUIUI 010! “I UiUIU'IUIUIU'IUImUIUIUIUI-thNN NNNNNNNNNNNNNNNN PH H mine-binuiucbciute AJNEOAJNOOADPtdF‘Ht‘F‘Ht‘ NH QmUl-FUNH UMP Finite Elements Theoretical Simulation of Self-heating/ Spontaneous Ignition Experimental Introduction Experimental Parameter Range Oven Design Determination of Surface Heat Transfer Coefficient Powder Density Particle Density Bulk Density Sample Preparation Time/Temperature Profiles The Differential Scanning Calorimeter DSC Equipment DSC Cell Temperature Control The DT Signal The DSC Signal The Mettler Processor Printer/Plotter Calibration of the DSC Determination of Tlag Sample Preparation and Insertion Analysis of Substances Theory 3 Development of a Finite Element Model for Self-heating [Self-ignition Introduction Mathematical Description of Heat Transfer Problem Finite Element Solution Finite Element Solution for Inert Product Finite Element Formulation Finite Element Grid Element Conduction Matrix Element Capacitance Matrix Element Convection Matrix The Global Matrix Equation Solution of the Global Equation Verification of the Heat Transfer Model Comparable Analytical Model Sample Problem Temperature Profiles Solution for Product with Heat Generation Inclusion of the Heat Generation Term The Force Vector Element Matrix tdii 78 79 83 84 86 89 91 91 91 94 96 98 103 103 103 106 106 110 111 111 112 113 114 115 115 116 117 119 120 121 123 123 124 125 125 126 126 127 127 139 139 141 6. Results and Discussion 6.1 Determining Surface Heat Transfer Coefficient 6.1.1 MIT Data for Avonmore Skim Milk Powder 6.1.1.1 Accuracy of Oven Measurements 6.1.2 The Kinetics of Self-heating 5.1.2.1 Density 6.1.2.2 Thermal Conductivity 6.2 DSC Data 6.2.1 Calorimetric Data for Milk Powders 6.2.2 Calculation / Preparation for Computer Simulation 6.3 Simulation Results 6.3.1 Detailed Simulation Run 6.3.1.1 Specific Heat 6.3.1.2 Activation Energy and Heat of Combustion 6.3.1.3 Simulated Time/Temperature Profiles 6.4 Simulation / Oven Validation 6.5 DSC / Oven / Simulation Comparison 6.5.1 DSC Accuracy 6.5.2 Conclusions 6.6 Sensitivity Analysis 6.6.1 Effect of Sphere Radius on Combustion / Ignition 6.6.1.1 Minimum Ignition Temperature 6.6.1.2 Time to Ignition 6.6.2 Effect of Activation Energy on Combustion / Ignition 6.6.3 Effect of Ambient Temperature on Combustion / Ignition 6.6.4 Effect of Surface Heat Transfer Coefficient on Combustion Ignition 7. Summary 8. Suggestions for Future Study AREEEDIQEE 1. Derivation of Element Matrices 2. DSC Programming Procedures 3. Computer Program Listings 4. Bibliography 144 144 148 152 153 156 158 158 159 161 163 164 165 165 166 168 175 178 179 181 182 182 185 188 191 194 199 202 206 215 218 235 ....... . N HUNHNHN H 01 UI Uih-bbUUH H mma‘ as as ea e e e 0 ”NH“ PH ‘0” ~IO\ LII-h O p mmmmm 001 mm mm . P'Hifl Sign: 6.15 6.16 A3.1 A3.2 A3.3 A3.4 LIST OF TABLES Typical case history of fires in a milk powder drying plant. Some major reported milk powder fire events. MIT data for milk powder cubes. MIT values for milk powder cubes. Powder sample radii and associated 'h' values. Time/temperature data for 4” sphere. DSC Calibration data. Analytical and numerical solution for the aluminium sphere. Analytical and numerical solution for a cooling problem. Time/temperature data for varying convection conditions. . Time/temperature and 'h' value data for 4' sphere. Experimentally determined values of the surface heat transfer coefficient. Radius/experimental MIT data for Avonmore skim milk. Bowes plot data for Avonmore skim milk. Experimental and derived milk-powder sample data. Calorimetric results for skim milk powders. Kinetic data for milk powder, derived by traditional and DSC techniques. Oven and model MIT values compared. Comparison of MIT ('C) values for different size spheres using different techniques. Simulation parameter range. 'Standard' simulation data. MIT vs sphere radius. TtI vs sphere radius. TtI vs Activation Energy. Simulated ambient temperature vs TtI. Heat transfer coefficient vs TtI. Programme MNFN3A. Subroutine SBFN3A. Subroutine FAT. Subroutine FBT. 58 60 96 99 113 128 133 136 145 147 152 154 158 160 162 169 178 181 182 182 187 190 191 196 218 219 225 226 A3.5 A3.6 A3.7 A3.8 A3.9 Subroutine FCT. Subroutine ARRAY. Subroutine SIMQ. Subroutine DQG32. Subroutine GMADD. A3.lo Subroutine GMPRD. 227 228 229 232 233 234 Hid . O NH U new UN H hhbfihéhbuu .......... @QGU‘bUND—‘U‘fi U! mmmmmmmm @QOiUI-huNH m 01 LIST OF FIGURES Schematic Diagram of Typical Drying Plant. 8 Heat Generation vs Meat Loss, Stable and 13 Runaway Reaction Zones. Determining Kinetic Parameters with F-K Parameters. 34 Heat Transfer Models for Self-ignition. 35 Typical Conical Spray Drier Section Showing Temperature Range and Powder Deposit Areas. DSC Enthalpy Range. 59 DTA and DSC Compared. 71 Test Oven. 75 Aluminium Sphere Heat-up Curve. 88 Pycnometer. 90 Jolting Volumeter. 92 Milk Powder Ignition Curves. 93 Typical DSC Printout. 104 DSC Cell. 105 Temperature in the DSC Measuring Cell as a Function of Time. 105 T Signal (Primary Signal) as a Function of time, T, and the Reference Temperature Tr. 107 n is the Vector Normal to the Sphere Surface. 118 Finite Element Grid. 121 Analytic vs Finite Element Solution for Aluminium Sphere Heat-up. 130 Analytic vs Finite Element Solution for Aluminium Sphere Heat-up (Expanded View). 131 Analytical vs Finite Element Solution for Apple Cooling. 135 F.E. Technique used to Model Different Values of Heat Transfer Coefficient. 138 4" Aluminium Sphere Heat-up Curve. 146 3” Sphere Milk Powder Ignition Curve. 149 2" Sphere Milk Powder Ignition Curve. 150 1.5” Sphere Milk Powder Ignition Curve. 151 Bowes Plot 155 Simulated Temperature Profile for 4" Sphere. 167 Predicting MIT Using Oven Data. 168 Oven and Model Profile at Lowest Oven 171 Ignition Temperature (1.5" sphere). 6.9 6.10 6.11 6.12 6.13 6.14 6.15 6.16 Oven and Model Profile at Lowest Oven Ignition Temperature (4” sphere). Oven and Model Profile at Lowest Oven Ignition Temperature (3" sphere). Oven and Model Profile at Lowest Oven Ignition Temperature (2" sphere). Predicting MIT Using DSC Data. Simulated MIT as Radius for Milk Powder. Simulated Time to Ignition vs Radius. Simulated Time to Ignition vs Activation Energy. Simulated Time to Ignition vs Ambient Temperature. Simulated Time to Ignition vs 'h'-value. )dii 171 173 173 177 184 186 189 192 195 1W Self-heating and spontaneous ignition are phenomena usually associated with a number of different situations. Firstly, in the case of self-ignition of hay-stacks, it may be linked to the exothermic biodegradation reaction occurring in the hay (Rothbaum, 1963 ). Alternatively, spontaneous ignition leading to a dust explosion can be seen as a natural occurrence in the coal industry, due to the dusty, inflammable nature of that product. What have gone largely uncatalogued until recent times are the many cases of self-ignition which occur in the food industry. In the food sector, the two conditions already referred to (i.e. exothermic reaction. and/or inflammable jproduct ) exist widely, though perhaps not on the same scale as the large heat generation due to the biodegradation in a damp hay stack or the highly inflammable dust clouds in a coal handling operation. Many foods are produced today in powder form , or have powder as a byproduct or as a waste material. In many cases, the chemical interaction of the powders with oxygen is a low-grade exothermic reaction. Unless the heat generation is dissipated to the environment, the product temperature will steadily increase over time, causing the heat generation term itself to increase . Either by the mechanism of internal heat generation or through the external supply of heat (e.g. a hot surface in contact with the powder or electric / friction spark) the food system may proceed to a runaway thermal reaction, i.e. a fire and / or a dust explosion. Instances of the 'food-based' spontaneous ignition incidents have occurred in such diverse areas as the dairy industry, with a wide variety of milk-based powders being recognised as potential fire hazards, and the confectionery trade, where sugar / chocolate based formulations pose a similar hazard. The cereal / grain business must also contend with fires caused by this phenomenon. Some fatalities have resulted from fires in the transportation and handling of grains in siloes, and in the unloading of grain cargoes from ships. The aim of the present study is to determine the conditions which give rise to potentially serious self- heating in milk powders (i.e. self-heating of powders which, either in storage or in the production line, proceed to a self-ignition situation). The International Dairy Federation (1987) summarised the danger by stating that W 1W. Self-heating / self-ignition in dairy powders is a major source of concern to the dairy industry at large. In -2- Ireland alone, there is an average of one major fire reported each year. Due to the competitive nature of the industry there is no doubt that there are a number of other fires or 'near misses' which go unreported. Typical damages from such an incident exceed $3 million, when equipment and building damage are totalled together with lost production capacity and lost markets (Wyeth, 1980). Serious injury to personnel in the dairy industry has fortunately been not been a facet of any of the reported fires to date. Fatalities have, however, happened in other branches of the food industry as a result of this hazardous phenomenon, adding greatly to the final toll of damages. Along with the major fires there are also numerous other smaller fires or occurrences of 'burnt particles' in spray dryer plants. The figures quoted here for Ireland are also . typical for the United Kingdom, and for Europe in general. One company , a market leader in detection and suppression of the explosions which typically follow a powder fire , reported a total of over four thousand dust explosions in Europe in the last twelve years, quite a number of which were in. the dairy sector (Graviner, 1990). Bartknecht (1989) observed that there is on average one industrial explosion for each working day in the industrialised countries of western Europe. Milk powder fires are difficult to predict or prevent. The prediction techniques currently in use are discussed in detail below, along with their inadequacies. Conventional fire detection is based on either flame detection ( using an ionisation chamber ) or smoke detection ( using a photocell ). In the inherently dusty environment of a spray dryer or conveyor tunnel, in which strong air currents are part of the drying operation, such devices are ineffectual. They are prone to much 'nuisance tripping', thus causing unnecessary equipment downtime and production losses, which most production schedules can not tolerate. Confidence is poor in such systems and production personnel view them as an obstacle to efficient performance rather than a vital production aid. Corporate policy in some companies thus prefers to suffer the occasional loss of a dryer installation due a fire ' rather than to absorb the recurring expense of driers shutting down periodicalry as a result of over-sensitive or faulty detection devices. For safety reasons, these devices err on the side of safety, i.e.,they must be more rather than less sensitive. The other detection parameter available is the dryer ambient temperature. As this temperature is used as an operational parameter it is not effective to combine this role with that of fire detection parameter.Thus an increasing dryer temperature causes the dryer control -4- system initially to increase the rate of product input. A critical threshold must be delineated to assess when this increased temperature is the result of a fire, rather than a routine operational deviation in the heat/mass balance. The main sensing parameter used in spray dryers , in the context of fire detection, is the change in chamber pressure when an explosion occurs. This is used to trigger a set of vents which are designed to funnel out the explosive pressure from the chamber to designated areas outside the dryer, where damage to buildings and/or danger to personnel are minimised (Bartknecht, 1989). Fire/explosion suppression may also be triggered by pressure sensing using either inert gas or steam flushing of the chamber to smother any incipient fires. Other techniques with potential applications include the detection of carbon monoxide, the presence of foul smell, infra red detection and a regular scorched particle test (IDF, 1987 ). The aim of the present study is to prevent the powder entering the fire or explosion phase. A number of reasons have been identified as the 'immediate' causes of milk powder fires. These include : (l) Self-heating of deposits, (2) External or friction heating, (3) Equipment malfunction, and (4) Start-up ~ conditions : damp dryer/product, improper heat/mass -5- balance (Synnott et al.,1986) The common denominator among these potential causes is that they involve some form of either excessive heating (from an internal or external source ) or restricted heat dissipation. This raises the temperature of the product to a point where the 'low level' exothermic oxidation reaction proceeds rapidly to a fire situation. IEDHEIBI_£QEI§£1 A typical case history of self-heating of milk powder has been described by Beever ( 1984 ). In a U.K. dairy plant manufacturing a range of milk powders, eight fires were reported in two and a half years. A schematic diagram of the plant is shown in Figure 1.1. Table 1.1 summarises the details of the fire incidents. Damp product features as an important contributory factor in the instances cited here. Start-up conditions play an important role as well since at start-up an imbalance usually exists between liquid / atomised product entering and the amount of heat being supplied to the dryer. This leads to the occurrence, in practically all commercial start-up situations, of burnt powder particles in the first run of product. These burnt particles , depending on their size and temperature, may, in conjunction with other process / environmental conditions, form the nucleus of a fire which may threaten -6- the entire powder area from the dryer through to the storage siloes. T§h12_111 Typical case history of fires in a milk powder drying plant, (Beever, 1984). DATE PRODUCT INCIDENT 8/20/80 22% Fat filled White spheres, (up to 5 cm in diameter) charred inside, found in sieve. 8/27/80 26% Fat filled Explosion in dryer. 12/20/80 26% Fat filled Burning smell; glowing patches on walls of dryer. 4/25/81 Skim milk Burning smell: burning spheres in sieve. 8/21/81 22% Fat filled Fire downstream. 10/6/81 26% Fat filled Charred lumps downstream. 2/23/82 32% Fat filled Fire downstream. 2/24/82 Skim milk Fire. It is difficult to extract accurate information on commercial milk powder fires as companies tend to maintain a veil of silence around such incidents. Fires may also damage a lot of the evidence although a case history such as that listed above generally indicates a history of 'near misses' occurring prior to a major disaster. Some fire insurers now insist on a full record being maintained of the temperature profile in a dryer, particularly if a fire / explosion has occurred with a -7- To Atmosphn Cancun“ r——’—i 0 WE l Cyclone i Air 811111 To Chlhr A my Dn'cr cm m "' >’_ I l Exhaust t. PM U Pins Exhaust Air nu rim V Instanti'zor No .2. Fins Collection Prior to Return to Chamber Bagging Fig. Ll : Schematic Diagram of Typical Drying Plant. '8- particular product / technology, (Golden Vale, 1984). In a survey of milk powder plants in Ireland twelve incidents of fire or explosion were reported over a six- year period, ( O'Callaghan et a1., 1988 ). Table 1.2 lists some of the major fire incidents reported in the Irish dairy powder sector in the eighties. Igblg_1;z Some major reported milk powder fire events. 2392291 IHQIQEEI Infant Marriott-Walker Box Drier chamber destroyed: Formula Electrical spark suspected. Smouldering in storage silo : Burnt powder clump as heat source of (WYETH). Skim-milk Major, slow burning silo fire (AVONMORE). Skim-milk Niro dryer fire (NCF). Fires occur both in the conventional conical type of dryer as well as the Box-type dryer. The Box-dryer has the additional risk of the filter bag area where very fine powder accumulates and where, for personnel safety reasons, electrical alarm buttons must be located . Storage siloes are also likely to present fire risks to the producer (Avonmore, 1987). While Table 1.1 does implicate the wet cleaning procedure ( where fires start occur on start-up after a CIP cycle ) , companies who operate a dry-clean procedure (e.g. Wyeth) have also had a -9- catalogue of fires and 'near misses' for other reasons (Wyeth, 1982). 1;1_IEE_QQBI_IQ_IED!§IBI The impetus for the present work comes from a strong demand from the Irish milk powder manufacturers to obtain a better understanding of the causes of fires in spray dryers, conveying systems and powder storage siloes. The phenomenon of self-heating in powders, and the related problem of dust explosions, constitutes a major hazard in the dairy industry. This has major implications for Ireland where dairying is one of the key indigenous industries. Dairy exports in 1987 were valued at I£1.2 X 109 ($1.9 X 109), produced by an industry employing almost eight thousand people in the manufacturing sector. Over 36% of these exports were in powder, ( An Bord Bainne, ,1988 ). Much of this product is manufactured in large central processing facilities which are similar to modern plants in the rest of Europe. With an increasing tendency to amalgamations, there is a greater dependence on central ised , large-volume facilities. These facilities are geared for longer running times and more automated operation. Hence the risk is increasing and there is a greater urgency to solve the problem. The European Community's intervention policy -10- which gave rise to the so-called 'milk-powder mountain' means more powder in storage, in bigger siloes, leading to greater risks of self-heating in the powder stocks. In Ireland alone, there is, on average, one major 'powder fire' reported every year. Considering that Ireland currently produces 5.2% of the European Community's milk quota , the investigative work in the present project has wide application in the dairy industry across Europe and worldwide. At current costs, a fire in a spray-drying chamber amounts typically to $3 million in damages and lost production (Wyeth,1980). The frequency of such incidents is similar in the.United Kingdom and in mainland Europe (Institution of Chemical Engineers, 1977 ). In the United States there is a substantial dairy sector, 'producing over $45 x 109 worth of dairy products per annum. Spray drying makes a major contribution in this area with over one thousand plants manufacturing dairy powders (U.S. Dept. of Commerce, 1990). In Michigan, for example, some 59 x 106 tons of non-fat dry milk are produced annually, approximately 6% of the national production (USDA, 1989). The problem of self-ignition of milk-powder is thus expected to be a significant one, both in the State of Michigan and in the 0.5. as a whole. ~11- 112.1333HQDIEAEIQE A brief summary of the heat balance involved in this form of thermal runaway reaction is shown in Figure 1.2. Two heat terms are involved -G : (1) Heat generation within a sample clump of powder, and (2) L : Heat lost from the surface of the powder. If G < L, the reaction is stable, since the generated heat is capable of being dissipated to the environment. If, however, G > L, then the heat will build up within the powder and. the reaction. proceeds from self-heating to self—ignition. The work reported in this thesis attempts to predict the threshold conditions in the balance of heat generated vs. heat lost. Commercially manufactured powders are used and the test and model conditions reflect, as far as possible, the conditions found in the commercial environments being modelled. The techniques currently used to identify the commercial conditions of self-heating are laboratory-based and empirical in nature. The model developed here uses a specific calculating procedure to follow the temperature profile in the powder over time until the point of thermal runaway is identified. With the programme it should be possible, with knowledge of the physical/thermal properties of a particular powder, to accurately predict -12- Thermal Runaway I. -> Rate of Heat Loaa G -> Rate of Heat Generation Temperature Fig. 1.2 Heat Generation vs Heat Loss Stable and Runaway Reaction Zones. c.13- the conditions to be avoided. This constitutes a significant advance over the existing experimental procedure which uses a 'bracketting' estimation of the Minimum Ignition Temperature (MIT); it is based on a series of laboratory tests, and thus has a basic in-built inaccuracy (Synnott et al., 1986). The MIT is classically taken as the mean of two temperature ranges . Using a suitable database of food properties, knowledge of product composition will be enough to allow the predictive model to operate successfully. While the model is developed for a spherical particle or a 'clump' of powder, extension to other geometries ( e.g. layers and cylinders ) is possible. The main thrust of the study is initially placed on the thermal conditions and the properties / composition of the product. The contribution of these factors to the likely onset of a thermal runaway reaction is assessed using both experimental data and the simulation. Several theoretical procedures are presently used in the literature to model self-heating of powders. To make the mathematics more amenable to solution, certain simplifying assumptions are made in the traditional techniques. These include the assumption of lumped thermal properties and high thermal conductivity in the so-called Semenov model. The assumption of negligible surface thermal resistance and zero reactant consumption is incorporated in the Frank-Kamenetskii model. The program developed in this -14- study incorporates in the model all available data on the product. and. its environment. Thus, finite surface and‘ internal resistance to heat transfer are allowed for in the model. Provision is made for inclusion of thermal / physical properties and kinetic data of the product and the particular exothermic reaction taking place. The study of self-heating / self-ignition problems has developed along two streams. Firstly, there is the empirical approach where pilot-scale, laboratory-based studies have been used to estimate the likelihood of various materials igniting under conditions of thermal stress. This work has mainly been conducted on materials associated with the building / construction industry. Public safety and good engineering practice dictate that only materials which are slow to ignite are certified as acceptable as building materials. Thus, in the event of a fire occurring, the structure itself will be able to withstand the fire temperatures for a certain specified time before igniting. Foodstuffs have not been subjected to the same detailed classification as construction materials in terms of combustibility. The phenomenon only~ manifested itself in -15- the food area with the relatively recent development of large scale central processing and storage facilities. The problems associated with the storage and shipment of grains fall into the category of dust explosions ; fatal incidents involving grain handling continue to be a hazard, and sometimes a source of fatalities in this industry. In recent times, the occurrence of powder fires and explosions, with the attendant danger to both personnel and plant, has caused the insurers of plants and the statutory bodies concerned with safety in the work place to inntroduce new regulations. These seek to minimise the risk of fires occurring and, if they do occur, to reduce the danger to personnel ( Irish Dept. of Labour, 1987 ). The present work is a contribution to the prevention of these fires. Studies in fire suppression techniques are also on-going to deal with the outbreaks which may still occur in spite of predictive work or detection systems. Developments in explosion venting is a further aspect of damage-limitation research at the post- ignition end of this phenomenon (O'Callaghan et al., 1988). The development in the study of this thesis has been based on the two models previously referred to (i.e. the Frank-Kamenetskii and Semenov models). These have given rise to a number of 'standard' parameters used to estimate the possibility of a product / material entering on a -16— thermal runaway reaction. Two such parameters are : ( 1) W and (2) the 2122115192122: W. The definitions of these two parameters give an indication of the nature of the scientific-cum- empirical approach normally used in studying combustibility of powders. 141W Based on oven or hot-plate tests the Minimum Ignition Temperature (MIT) is 'defined' as the average. of the lowest ambient temperature at which a particular size and shape of sample ignites, and the highest at which it fails to ignite (Synnott et al.,1986 ). Bowes and Townsend ( 1962 ) suggest that the "high" temperature should produce at least five non-ignitions to be accepted as the upper bound. Thus, the MIT is weighted towards defining a . safe, conservative maximum operating temperature. This is in keeping with good engineering practise. The Frank-Kamenetskii parameter, Equation [1.1], has been described as a 'dimensionless reaction rate' and incorporates the thermal and physical parameters pertaining to the oxidation / combustion reaction under study [Thomas,1960 ] . Thus, -17- 6 a QfEArozexp(-EA/RT3) [1.1] kRTaZ where Q is the heat of reaction per unit volume, f is a frequency factor, EA is the activation energy, ro is the critical sample dimension ( radius of a sphere, half width of slab, half-side of a cube), R is the universal gas constant, Ta is the ambient temperature, and k is the thermal conductivity. 6 represents a threshold condition for ignition. Computation and interpretation is clearly dependent on many factors. Much of the theoretical work on this self-heating predates the solution techniques available on digital computers. It has been presented in a form which is not easily available ‘ or useful to those faced with design or operational decisions in ensuring the safety aspects of spray drying. Accordingly, it is well outside the scope of the those charged. with directing' the safe operation of a spray drying plant under varying environmental conditions and product specifications. The theoretical/empirical body of work recorded is not amenable to predicting possible combustion risks of the many food powders which are dehydrated today in spray drying facilities. Thus, the available literature is either highly empirical -13- in nature or is circumscribed with theory aimed at legitimising empirical / pilot—scale results. As a result of the present work, a 'user friendly' program is available which allows firm prediction of the combustion indices such as MIT and 'Time to Ignition' for a variety of dairy-based powders. This will assist the food engineering designer as well as the product formulator; and it should result in a major improvement in safety in the workplace. Simultaneously, this should mean the elimination of a source of major downtime and equipment losses. 21.931391IIIE The primary objective of this work is to predict the conditions under which milk powders undergo self-heating which proceeds to a thermal runaway reaction, either in the drying plant, conveying system or in subsequent storage in siloes or other bulk containers. These conditions tend to be particular to certain milk powder products or groups of products (e.g. fat-filled or non— fat-filled powders). Stated concisely, the study is directed to achieve the following goals : (I) Develop a mathematical model to predict the temperature profile in a spherical product sample with: a) Beat generation: b) A Newtonian boundary condition typical of a spray dryer: c) Variable thermal properties. (II) Develop reaction kinetics' and product property data to model the combustion / heating term using the Frank-Kamenetskii and DSC methods. -20- (III) Verify the self-heating/runaway reaction model experimentally. (IV) Identify the critical combustion parameters or parameter ranges for commercially-produced milk powders while in spray dryers and subsequent storage. -21— 2&-LIIIBLI!B:_B£!IE! The schematic diagram of a body undergoing an Arrhenius exothermic reaction and simultaneous heat loss from its surface is shown in Figure 1.2. Essentially if the convective heat loss, L,is greater than the heat generation term, G, then the system will cool down. If the reverse holds, i.e. G > L, then the reaction may proceed to an unstable runaway state (Drysdale,1985). 2I1_IBADIIIQEAL_IIE£_£_£QH!Q§IIQE_I§EQBI The bulk of the research on the phenomena of self-heating and self-ignition , experimental and theoretical, is based on two theoretical developments, those of Semenov (1928) and Frank-Kamenetskii (1939). Both theories have limitations arising from simplifying assumptions used to facilitate the mathematical solutions. The assumptions may or may not hold depending on the test conditions or products being tested. The theories and associated assumptions are outlined below. 11111.133_§£H£EQ!_HQD£L Semenov (1928) sought to predict a critical ambient temperature, Ta,cro above which thermal runaway occurs and below which the system eventually cools down. He equated -22- the heat generation term, G, with the heat being lost at the sample surface, L : a II I." [3.1] The temperature is assumed to be uniform at all times throughout the reacting sample. Heat lost at the surface depends on the convective heat transfer coefficient. The heat generation term G is according to the general Arrhenius reaction: a - Q cin r e’EA/RT [3.2] where the terms on the right hand side except C1“ are as previously defined for Equation [1.1]. C1 is the concentration of reactant component i , and n the order of the reaction. IL is equal to the heat loss from the surface and may be written as: L = h A ( T - Tamb ) [3.3] where h is the surface heat transfer coefficient, A the surface area, T and Tamb the surface and ambient temperatures, respectively. At the critical or threshold point shown in Figure 1.2, -23- c = L [3.4] Thus Q Ci“ f e‘EA/RT = h A ( T - Tamb ) [3.5] As L is tangential to the G curve, d9 - dL [3-6] dT dT at the critical temperature. Hence Ea Q cin f e‘EA/RT = h A [3.7] R T2 Dividing [3.5] by [3.7] gives R T2 / EA = T - Tamb [3.8] where now T = Tcr . Thus: R Tcr2 / EA = Tcr ’ Tamb [3-9] Rewriting: R Tcrz / EA - Tcr + Tamb - 0 [3.10] Using the standard solution for quadratic equations the threshold temperature Tcr has the following values : Tcr - { 1 + ( 1 - 4 T R / EA )1/2 } 2? A = EA/ 2 R { 1 + ( 1 - 4 Tamb R / E )1/2 } [3.11] Equation [3.11] specifies the product temperature above which a fire situation may occur. This condition exists, i.e. Equation [3.11] has a 'real' solution provided EA > 4 Tamb R. A typical milk powder has values of EA = 40 kJ / mole and R =- 8.34313 J / K mole. This condition gives a value of Tcr < 1199 K, a temperature clearly not exceeded during powder manufacture and storage. Thus ignition is theoretically possible. To further evaluate Tcrr the -24- square root term in Equation [3.11] is calculated using a Binomial series expansion ( Abramowitz et al., 1972): (l-4TambR/EA)1/2 g l-4TambR - 16Tamb2R2 - 64Tamb3 R3 2 EA 8 2A2 16 3A3 [3.12] From Equation [3.11] Tcr a EA/ZR { 1 - ( 1 - 4 Tamb R / EA )1/2} = Tamb + Tambz R + 2 Tamb3 R2 + ... EA 2A2 [3.13] In this technique, it is usual to evaluate the series expression in [3.13] by using only two terms. HenCe the critical temperature rise above ambient is: ATcr = Tcr ' Tamb = Tamb2 R /EA [3-14] The error involved in truncating this series has been estimated by Simchen (1964) for various combinations of temperature and the Activation Energy. With Tamb = 500 K and EA - 168 kJ / mole, the truncation error is equal to 5.0 %. Semenov's approximation assumed that -25- Tcr = Tamb [3.15] arcr = R Tcrz / EA [3.16] Equation [3.15] is the more correct version of the development based on the Semenov theory (Gray et al.,1967A). Equation [3.14] thus gives the maximum spontaneous temperature rise that may occur within the system without ignition of the sample. The principal limitation of the Semenov theory is the assumption of negligible internal resistance to heat transfer which is not valid since the thermal conductivity of milk powder is in the range 0.04-0.1 W / mK, (MacCarthy, 1983 ). The Semenov model is best suited to situations where good convective heat transfer occurs in a sample undergoing heating. Merzhanov et a1. ( 1961) found the Semenov theory to suit well the case of heating of a well stirred liquid explosive. In summary . Manes—W11 Mise- 111aZ_IEI_IBAEI:BAHEEII§EII_HQDIL While the Semenov model assumes that convective heat transfer at the surface is the only' barrier to heat -26- transfer, the Frank-Kamenetskii (1939) theory assumes that conductive heat transfer through the sample is the slowest, and hence is the dominant mode of heat transfer. This theory gives rise to a number of important parameters used in categorising combustibility. An unlimited supply of reactant is assumed for the heat generating reaction. Frank-Kamenetskii (F-K) assumed an exponential approximation for the temperature dependence of the reaction rate, kT: kT - A e’EA/RT [3.17] Here A is a lumped pre-exponential term. At the point of criticality, which is the area of interest , Equation [3.14] allows the exponent of Equation [3.17] to be rewritten as : EA / RT = EA / [R ( Tamb + 11)} - EA / [R Tamb ( 1 + AT/Tamb )1 [EA/R Tambll 1'AT/Tamb+AT2/Tamb2] [3.18] The last step is accomplished using a Taylor series expansion of the term 1/(1 + T/Tamb) ( Kreyszig, 1967 ). Near criticality, T << Tambz and hence only the first two terms of [3.18] are significant. Equation [3.17] now takes the form: -27- R Tamb R Tambz [3.19] Thus RT = kTamb e‘Ea/RTame [3.201 F-K solved the equation for heat conduction with heat generation de12 T + c = 1 0T k a at [3.21] For uniform symmetrical heating, this reduces , for the unidimensional case, to : 1221. + r 21 + c = .1. 21 or2 r or k a at [3.22] K is constant depending on the geometry of the sample under investigation. It has values of 0 for an infinite slab ( of thickness Zro ), 1 for an infinite cylinder (radius ro ) or 2 for a sphere (radius to) (Kreyszig,l967). The F-K model assumes at the boundary that the surface -28- temperature equals the ambient temperature. Thus , for a sphere: T = Tamb at r = ro, sphere radius, t>0 [3.23] Symmetry around the sphere center yields the second condition: dT/dr = 0 at r = o , t>0 [3.24] To facilitate the solution process, the temperature and distance are non-dimensionalised using 6 and z: e = EA(T - Tambl/R Tambz [3.251 and z - r/ro [3.26] Equation [3.22] may thus be rewritten in dimensionless form as : k RTame 026 + K 69 e -o cin f exp -EA 6 r02 EA oz2 2 oz RTamb 1-e*e [3.27] -29- where e* = RTamb/EA- If e*<<1, [3.27] may be further approximated by : 029 +.§ 23 = rozEAchin x exp(-EA/RTamb)exp(6) 022 2 oz kRTamb2 [3.28] or del2 e = - 6 exp(8) [3.29] where 6 = r02 EA o f cin e'EA/RTamb k R Tambz [3.30] Frank-Kamenetskii assumed high Biot number conditions , i.e. Bi > 10 [ Bi = hl/k ] whereby the heat transfer is limited only by the internal heat transfer rate (i.e. the thermal conductivity of 'the sample). Thus ‘tne__§uzfiggg -- s .1 z - tea . é ' .-~umeq - 9' 9" '. . . In a spray dryer producing skim milk powder, for example, the conditions give the following parameter values h : 14 W / m2, k :0.1 W / m K and r0 : 0.05 m, resulting in a Di - 7. Thus this fundamental assumption regarding the heat transfer conditions is not satisfied by the industrial conditions of interest in the present study. While this is so, the F-K theory is the basis of much of -30- the work done to date in self-ignition studies. The application of the F-K model is thus further outlined below. 1) The reaction rate can be described by a single Arrhenius expression, Equation [3.2]. 2) There is no reactant consumption: combustion/ degradation reaction is zero order in the reactant concentration. The heat generation is not limited by the amount of powder or oxygen concentration,i.e. Ci“ - 1 or n - 0 in Equation [3.2]. 3) The Biot number is large ( i.e. minimal surface resistance to heat transfer, h - infinity ). 4) The thermal properties are constant ( e.g. thermal conductivity, density and specific heat are independent of temperature ). The method essentially involves solving the unidimensional -31- equation for conduction heat transfer with a heat generation term making the basic assumptions that : [1] Surface temperature = ambient temperature. [2] Temperature is symmetrical about center. Putting aT/ct 0 in IEquation [3.27] signifies ‘that steady state conditions have been reached i.e. (heat generation - heat loss). This gives a temperature distribution ( see Equation [3.29] ) where 5cr: as defined in Equation [3.30], is dependent on values of EA, R, Q, f, C1“ and Ta , parameters which describe the 'oxidation' / heating reaction of interest . 'k' is the sample thermal property for the assumed mode of heat transfer ( thermal conduction in this model ) and r0 is the characteristic dimension of the system. Values of 6 greater than the critical value, 5crr will cause ignition and a runaway reaction. The temperature, 9, corresponding to the critical value 5cr: is the MIT value. 6 contains all the information about the reaction and the reaction conditions, and can thus be used to investigate how the various parameters, independently or in combination, influence each other. This is of particular interest when studying a combination of factors that cause 5cr to occur since this is the ignition threshold circumstance. 5cr is also useful to work with because it can be identified experimentally by the widely used oven and hot-plate -32- techniques [ Synnott et al., 1986 ]. Rearranging the definition of 6cr in Equation [3.30] yields : 0Ta,cr2 - EAQfCine(EA/RTa,cr) r20 kR [3.31] Taking natural logarithms : 1n(5cr T2a,cr/r20) = 1n(EAQani/kR)’ EA R Ta,cr [3.32] Defining constants M and N allows Equation [3.32] to be written as : 1n 5cr T2a cr = M - n [3.33] r20 Ta,cr where M a N are the intercept and slope, respectively, of the plot of ln (6chZa'cr/r20 ) against 1/Ta,cr as shown in Figure 3.1. -33- 1n (5cr2 Ta,cr2/roz) 1 / Ta,cr zigg:g_1&1 Determining Kinetic Parameters with F-K Model. H is defined as In ( EA 9! C13 / k R ) (Intercept) and u is a; / R (81096) The mathematical / graphical technique forms the basis of most of the research reported to date on cataloguing the self-ignition properties of powders. Figure 3.2A summarises the temperature profile assumptions inherent in the Semenov and F-K models. In contrast to these models, the temperature profile shown in Figure 3.2B includes finite surface and internal thermal resistances. The U.C.C. profile in Figure 3.2B is employed for the model developed in this study. It is more realistic than the other profiles as it simulates the actual situation during the spraydrying and storage of dairy powders. The -34- Surface T. Amb __ _— _T. Amb. SEHENOV FRANK - KAr‘iENETSKll 3.2 A T. AmD.——. .__T. Amb‘ U.C.C 3.2 8 Figure 3.2 Heat Transfer Models for Self-ignition. respective Semenov and F-K assumptions have been shown to be inapplicable to dairy powder ignition. Gray et a1. (1958) studied the applicability of the reaction rate approximations to their theoretical solutions for the Semenov and F-K cases. They concluded that two equations give good approximations to exp(- EA/RTamb) in estimating Tcr for temperatures within the range 0 < e < 2, in the extreme F-K and Semenov cases. The first; was the Frank-Kamenetskii type "exponential" approximation whereby expt-EA/RTamb )= exp(-EA/RTamb)exp(9) [3.34] where G is a dimensionless temperature ( defined in Equation [3.25]). A quadratic approximation is also possible i.e. exp(-EA/RTamb) = {exp(-EA/RTamb))(1+0.726+62} [3.35] Gray et al. also considered an approximation where the heat transfer equation is solved using an "average" cm' -35- lumped internal temperature value, as suggested by Rice et al. (1935). To determine the ignition-induction time Rice et al. assumed that the heat lost during an explosion is negligible compared. to the generation of heat in the chemical reaction. This method calculated induction times by equating heat generation with heat storage, ignoring both finite conduction and convection effects, and thus limiting the general applicability of the analysis. 5cr Equation [3.30] defines the Frank-Kamenetskii critical parameter, the 'dimensionless reaction rate', 6 or 5cr- It is used to assess the combustibilty of different products. For a given product type and sample shape / size it signals whether or not a sample will ignite. Beever (1984) observed that 5cr depends mainly on sample geometry and to a lesser extent on sample material. Gray et al. (1967) reviewed the different methods used to determine criticality and expressed 5cr as : 5cr -( Bi / exp ) ro ( s / v ) [3.361 where S and V are the sample surface area and volume -37- respectively. In summary: 5cr = Constant x Heat Exchange Factor x Shape Factor [3.37] Equations [3.36] and [3.37] summarise the variability of the 5cr factor, on which the theory and MIT estimates are based. The environmental conditions, product thermal properties and the physical size and shape of the sample influence the value of 5cr which is often taken as a constant. Thomas (1957) examined the effect of heat transfer conditions, i.e. theBiot number, on 631-, and graphically showed 5cr as a function of Bi. He proposed asymptotic 5cr values of 3.32 for a sphere, 2.00 for an infinite cylinder and 0.88 for a slab, as Bi becomes very large and - approaches infinity. This set of values is equal to those proposed by Frank-Kamenetskii (1939) for the use of the criticality factor in samples of standard shape undergoing symmetric heating / cooling. Walker (1961A) , using data derived from studies of spheres of pie wool, claimed that 5cr is temperature dependent and that these values of 3.32, 2.00 and 0.88 are lower than the values which occur under normal ambient temperatures. Thomas et al. (1961A) proposed ‘values for’ spherical specimens of ‘wood fibre insulation between 1.7 for smaller samples (ro < 1/8 ") -33- and 3.2 for larger samples (e.g. r0 > 4" ). Other values have been proposed for different products. For skim milk, typical values are : 3.46 for a sphere, 0.92 for a layer, 2.64 for a cube and 2.89 for a cylinder ( Beever 1984, Synnott et al. 1986 ). Mathematically’ derived 'values for 5cr for' the sphere, infinite slab and infinite cylinder have been reported by a number of workers (e.g. Boddington et al (1971)). The proposed values for 5cr were , respectively, 3.322, 0.878 and 2.00. Walker et al. (1975) discussed the effect of temperature on the 5cr values, plotting 5cr values as a function of EA / R Tamb‘,cr . They also recorded the inherent instability of 5cr values when the temperature rise prior to thermal explosion is relatively compared to the absolute value of ambient temperature. This prerequisite is built in to the development of the critical parameter by Frank-Kamenetskii. Truncating the Binomial series in Equation [3.14] implies that EA was temperature dependent. Using a numerical technique, Anderson et al. (1974) similarly arrived at values of 3.322 for a sphere and 2.0 for an infinite cylinder. The 5cr parameter is still a widely-used index of combustibility. When it is used in conjunction with Equation [3.31] and graphs such as Figure 3.1, extrapolation is possible between different shapes. -39- O'Mahony et al. (1986) investigated the effect of sample shape and size on self-ignition of a fat-filled powder. They found that extrapolation from results of spherical samples gave very good correlation with experimentally derived results for layers of material. These latter results were achieved using 'hot-plate' tests on milk powder. Extrapolating from results on cubes did not produce good correlation. This is an important result as industrial fires frequently begin when layers of powders form in the dryer, conveying system or silo ( Figure 1.1). Experimental and theoretical research conducted with samples of spherical geometry , as for the case in the work being reported here, thus has important industrial significance. Thomas (1958) expanded the Frank-Kamenetskii model to allow for the effect of Newtonian cooling at the surface. This boundary condition is similar to that represented in Figure 3.2B and is stated mathematically as : h ( T - Tamb ) = k ( dT/dr ) , r = r0 [3.38] Using the same development as outlined above to arrive at Equation [3.19], Thomas studied the effect of changes in the parameters h, k and r0 on the value of the Frank- -40- Kamenetskii critical parameter, 5cr~ He concluded. that where thermal resistance at the surface is significant, the standard values for 5cr (i.e. 0.88,2.0 and 3.3 for the slab, cylinder and sphere, respectively) may be too high. When the factor hro/k is small (i.e. Bi < 0.5 ) it is safe to assume uniform temperature within the material. This combines both simplifying assumptions into the one solution process but one of no practical relevance. A similar mathematical treatment was proposed by Chambre (1952). Thomas (1960) summarised the theoretical approximations of Frank-Kamenetskii and Semenov. He proposed a quasi- stationary solution for inert materials and an effective heat transfer coefficient Beff : Beff = (6cr X exp) / ( 1 + k ) [3.39] Values of 5cr were compared for standard shapes as obtained using three different approximations : 1) Exact solutions based on an exponential approximation to the Arrhenius term: 2) Approximations which account for conduction in inert material: and 3) Approximations using effective heat transfer coefficient term. He also considered the validity of the 5cr values defined using the effective heat transfer coefficient in the -41- context of possible reactant loss near the critical condition, and concluded that the values were applicable. The results obtained were in quite good agreement with those of Thomas (1958). Thomas et al. (1961B) applied the Frank-Kamenetskii theory to a reacting slab one side of which is in perfect thermal contact with a hot surface and the other side is exposed to a constant cooler temperature. These conditions closely simulate the environment in a spray dryer. In particular, it mirrors the situation where layers or clumps of powder occur on the hot dryer walls or floors. The standard approximation for the Arrhenius. reaction rate term, as used in Equation [3.18], was seen to introduce an error of less than 9% for a value of RTp/EA equal to 0.04, where Tp is the temperature of the hot surface. The critical parameter, 5cr: was presented as a function of hro/k and the ambient temperature. Thomas (1972) summarised the basic theoretical techniques to model the phenomena of self-heating and self-ignition. The Semenov and Frank-Kamenetskii models are reviewed: a dimensional analysis approach is used to arrive at the same results as these models. Ignition temperatures as derived either experimentally or theoretically are functions of the material properties, sample shape and size and the environmental conditions. Nominal or -42- effective values of the Activation Energies can be used in simple models to correlate the experimental results. Induction ‘times 'were found to depend on sample size. Although tests on product samples may show low or borderline values of heats of combustion, the low values can become significant under bulk storage conditions. Thomas (197) also studied, reviewed. and. described ‘the theory on self-heating and ignition. He also developed a number of approximations to extend the F-K theory, in particular to non-F-K situations. The work is based however on the two extreme models of Figure 3.2. The inherent assumptions are erroneous and do not hold when analysing the problem of self-heating of milk powders. The approach tends to be empirical, and is unsubstantiated by laboratory results. Boddington et al. (1971) proposed an approximate technique to determine criticality in bodies of arbitrary shape. The procedure is tedious for irregular geometries. They recommended equivalent sphere radii for shapes studied under 'Frank-Kamenetskii conditions' (i.e. high Biot number). For arbitrary Biot number, an empirical approach is needed. They found that expressing shapes as equivalent -43- spheres produces meaningful results. The technique is applicable to any shape with a center of symmetry such that the entire surface is visible from the center. The justification of the derivation of the criticality criterion is plausible but not rigorous. Gray et al. (1967B) also studied equivalent spheres and reported that, under Frank-Kamenetskii conditions, no body has a lower critical temperature (Tcr ) than a sphere of the same volume. Boddington et al. (1981) returned to their 1971 paper and outlined its application. Under the Semenov condition of low Bi number, they claimed ignition to be less likely as the reactant mass size is reduced. This is expected sinces with a small sample, the heat is dissipated faster from the center/ inner layers of product, minimising the danger of a hot-spot developing . A maximum error of 5% was reported for the models used . Where a small range of sample sizes is used, however, potentially dangerous extrapolation errors are possible in estimating criticality. The results of a finite difference solution of the non-stationary heat transfer problem in a sphere are also presented. Errors of less than 10% are claimed when the results were compared with similar steady-state solutions. Walker et al. (1975A) used published values of 5cr to -44- extrapolate ‘their results to other' shapes. This is a departure from the method previously outlined by Walker(196lB) where the same author advised against using the Frank-Kamenetskii parameter. The conditions modelled empirically by Walker et al. include infinite surface heat transfer and a straightforward zero-order reaction. The calculation of 6m. values depends on establishing 'equivalent spheres' for' the shapes for 'the thermal conditions prevailing just prior to ignition / criticality. As the heat transfer term is specific, the application of the model is limited. An extension of the 'arbitrary shape' model was reported by Walker et al. (19758) to account for variable surface conditions. Accuracies to within 3% were calculated for 5cr values, when values were compared to previously calculated values for standard shapes. The equivalent spheres used for the various shapes are those calculated for negligible surface resistance to heat transfer. Wake et al. (1976) produced results based on the use of a variational method by Wake et al. (1973) for which they claimed accuracies to within 0.1% for Class A geometries (i.e. infinite slab, infinite cylinder, sphere). The accuracy is with respect to the figures such as those of Thomas (1958), whose results were questioned by Walker et al. (19618). When compared with the results of a finite element treatment (Anderson et al., 1974), the accuracy at -45- ‘finite' Biot numbers falls to 24%. Confidence is shown by these authors only for conditions of infinite Biot number. No [experimental results were furnished. to support the model. An approximate method of estimating the influence of the ambient temperature on the criticality factor, 5cr: has been proposed by Walker et al. (1977) for all values of the Biot number, for regular and irregular shapes. For class A geometries, the method of Walker (1961A,B) was employed, without reference to the doubts the author casts on this method in the original publications. The results are close to those produced by the method of Simchen (1964) for values of EA / RTamb < 0.1. The work of Walker et al.(197SB) was used to extend the method to bodies of arbitrary shape. 11ZI1_IABIIE§_IE££!AL_£QEDEQIIZIIX Walker et al.(1978) referred to the theoretical development of Carrie et al. (1959) and the work of Boddington et al.(1971) to investigate the effect of varying thermal conductivity on the onset of criticality in reactions of zero order. An exponential temperature dependence was assumed for the thermal conductivity parameter. They concluded that both positive and negative temperature coefficients of thermal conductivity can give -45- rise to ignition. A negative exponential coefficient of thermal conductivity accelerates criticality: a positive coefficient retards it. Where the coefficient is greater than or equal to that of the chemical reaction rate (i.e. that of the heat generation term), thermal explosion does not occur. The same procedure ‘was followed by ‘Walker (1980) to assess the effect of an assumed linear dependence of thermal conductivity on temperature. He used a numerical technique to estimate 6 as a ratio of the conductivity coefficient and the exponential coefficient of chemical reactivity. The ratio defined Scr and #crr the dimensionless central temperature rise. The results are in good agreement with Walker et al. (1978) but both studies are only relevant for the Dirichlet or Frank-Kamenetskii boundary condition shown in Fig. 3.2A. MacCarthy (1983) proposes the Maxwell-Euchen model as the best suited for estimating the thermal conductivity of milk powders. An effective thermal conductivity, ke, is proposed : ke ‘ kair 1 ' fV( 1 “biksol/kair} ) 1+fv(b-1) where, kair a thermal conductivity of air fv - volume fraction of solid ksol - thermal conductivity of solid component -47- [3.40 b = 3kair/( 2kair + ksol ) Using a guarded hot-plate method, thermal conductivity values in the range 0.036-0.109 W/mK are reported for skim milk powder samples. This is the range of values accepted and used in the present model. 1e212_BIAQIAEI_QQ!§!HRIIQH The standard assumption used in ignition studies is that the reaction is not limited by availability of reactant present. In ‘the case of industrial scale self-heating reactions this assumption is reasonable since there . is always ample suply of both powder and oxygen available for the reaction to proceed. In practically all the theoretical studies previously conducted, the assumption of zero-order reaction is basic to the solution / analysis procedure. Some researchers have tried to include consideration of a non-zero order reaction, in particular with regard to analysing laboratory ignition tests, where the reduced scale of the tests implies that the reaction is not independent of concentration. Non-zero-order reactions result in additional complexity in the mathematical treatment. While depletion of reactant is significant after ignition and may cause the reaction to slow down or stop, the situation prior to criticality is approximately zero-order, even with the relatively small -48- samples used in laboratory tests. This zero-order assumption takes on a new dimension, however, when predicting MIT from tests on mg scale samples in the Differential Scanning Calorimeter. For the Frank-Kamenetskii boundary conditions, Tyler et al.(1965) calculated values of the 5cr parameter and times to ignition for the cylinder and sphere for zero, first and second order reactions for different values of temperature rise. Thus, for a given set of conditions,6cr has the standard value of 3.322 for a sphere, whereas it varies from 3.405 to 5.309 as the order of reaction increases to first and second order, respectively. Boddington et al. (1977) made a comparison between the SemenOV' and Frank-Kamenetskii treatments for self- ignition, with and without reactant consumption. They used a quadratic approximation to the Arrhenius term to obtain an approximate analytical solution. When reactant consumption was taken into account, the main difference in assessing criticality was an absence of a major jump or step function in temperatures at criticality. Essentially, the Boddington et al. paper is concerned with extreme conditions given the inherent strong simplifying assumptions made. Thus further work needs to be done to extend the theory to real industrial situations. The intermediate' / finite conditions shown as Figure 3.2B -49- above as the focus of the present work is an attempt to advance the theory in this direction. In summary, the major limitations on the Frank-Kamenetskii treatment include the restriction on the heat transfer model and the question of reactant consumption. Limiting the heat transfer to situations of high convection at the surface is an unrealistic condition, both for the test set-up in the experimental conditions described below and in many industrial situations. The role of reactant consumption is not as critical, as the bulk of the powder reacts during the ignition phase. Hence, the course of the reaction up to ignition is not effectively limited by the concentration of the reactant. Frank-Kamenetskii's assumption of a zero-order reaction may thus be taken as a good working assumption up to ignition, the phase of the reaction of interest in the present study. Powder fires or fires resulting from spontaneous ignition are not unique to the dairy industry. Studies of the phenomenon have been conducted in other materials as well. Products studied include wool, cotton, grains, insulation, coffee, iron filings , wood, sawdust, coal, fishmeal and tobacco. A brief summary of the various studies has been reviewed by Bishop (1981). -50- Thomas et al. (1961A) used the traditional self-heating self-ignition theory to analyze the results of Mitchell (1951) concerning ‘the thermal behaviour' of 'wood fibre insulating board. They concluded that the simple model, based on a single Arrhenius rate reaction, extrapolates well to other temperature and size ranges. They cautioned that both self-heating and self-ignition tests need to be conducted on a sample in order to derive the best overall picture from small scale tests. In studies on smouldering sawdust, Palmer (1957) found that there is a minimum layer depth below which smouldering' self-extinguishes. As the air 'velocity increases over a wood pile, the minimum depth decreases considerably due to the increased rate of burning at the surface. Particle size also affects the minimum depth, causing it to increase with increased particle size. Synnott et al. (1984) published data on the MIT values of self-raising wheat flour, with a moisture content of 11.7%, using the standard hot-plate method. Results varied from an MIT of 311‘C for 5 mm layers to 265'C for 20 mm and 238'C for 40 mm. The authors concluded that the "Bowes fitting equation", Equation [3.33], fits the results well. They also noted that best results are obtained for samples with a diameter-to-depth ratio greater than 2.5:1. Extrapolation to smaller layer depths gives rise to -51- problems since the smaller layers tend to buckle and crack under heat. Drysdale (1980) published a detailed review on smouldering combustion. Products discussed include sawdust, cellulose and polyurethane, polyisocyanurate and phenol formaldehyde foams. The transition from smouldering to flaming combustion is discussed. In haystacks, for example, the initial self-heating is due to the action of thermophilic bacteria. Chemical oxidation in a stack can raise the material to smouldering temperatures after which combustion may occur (Rothbaum, 1963) . Drysdale stated that for a material to smoulder, it must form a rigid char: such is the case with dairy powders. - '00 When wool is oiled by certain textile treating oils or contaminated by mutton tallow, an exothermic reaction may occur with atmospheric oxygen possibly leading to spontaneous ignition. This may occur even in dry sterile conditions due to the oxidation of the fat ( Walker et al., 1965A ). Walker (1961A) studied of the heat balance in spontaneous ignition at the point of criticality. He proposed solutions to the 'simpler' heat transfer problems as -52- depicted in Fig. 3.2A (i.e. either zero internal or surface thermal resistance), for regular geometries. He introduced a parameter 'y' which he defined as the rise in temperature needed to double the reaction rate under the conditions of the experiment, as an alternative to assessing the role of activation energy. The technique is based on Carrie et al. (1959); the results are independent of Arrhenius or other dependence of the reaction rate. Walker contradicted the commonly accepted view that 5cr is constant and his calculated values show 5cr to increase with temperature. Thus for a sphere, 5cr increased from 3.32 for RTamb/EA of 0.0 to 4.46 at RTamb/EA of 0.2. His technique avoids the use .of the approximation made in Equation [ 3.19 ] which is basic to the Frank-Kamenetskii approach and which the author claims enforced an unaccountable temperature dependence on both EA and kTambr Equation [3.20]. He allows for a variable thermal conductivity by using an 'average' value for the parameter. He also uses a temperature intermediate between central and ambient temperature as the 'average ' body temperature. A positive temperature coefficient of reaction was needed for spontaneous ignition to occur. The greatest source of error in the theory, as in the Frank- Kamenetskii theory, is the assumption of a zero order reaction and the inaccuracy inherent in using the 'average' temperature. The author estimated the relationship between the critical size and the surface .53. temperature in piles of pie wool where the surface temperature is equal to the ambient temperature. Walker (19618) developed the theory to include the case of a variable surface heat transfer coefficient. He claimed that this technique is an improvement on the traditional developments of the Frank-Kamenetskii / Semenov theories as calculation of the complex 6 parameter is not required. Neither do the constants in the Arrhenius expression of the reaction rate need to be determined. The method instead depends on the experimental determination of the reaction rate and the nature of the reaction's temperature dependence. The model of Thomas (1958) is only valid near absolute zero temperature. The author did not assert great confidence in his own model either, however, since the surface heat transfer coefficient can not be accurately determined for the model and the model equations are not stable for variations in the coefficient value. The model is further limited to reactions of zero order. Walker et al. (1965A) reported results that show that wool is more likely to ignite than hay or fibre insulating board. They studied wool, packed into spherical flasks, placed in a constant temperature bath at temperatures between 90‘C and llO‘C, in an oxygen atmosphere, for times varying from 1 hr to 200 hrs. They proposed a simple heat generation equation -54- G = 2.4 x 10-6 (20.008T / t0.282) [3.41] where G is rate of heat generation (cal/sec g) , at temperature T ('C) after time t(sec) . The equation does not include the Activation Energy, unlike many of the Frank-Kamenetskii based theories. The equation was used to calculate the heat generation in cylinders of clean and greasy wool in a dry air atmosphere (Walker et al., 1967), at temperatures between 95‘C and 160'C, based on the center temperature rise. Experimental results showed centerpoint temperatures at ignition of 160°C for clean wool and between 143-147'C for greasy wool. Further experiments were conducted on similar wool samples . by Walker et al. (1968) with the cylinders suspended in a heated oven. When the wool samples were predried ( by air 8 80'C for 8 hours ), they were found to ignite a few degrees higher than similar samples dried in nitrogen in an oil bath and then exposed to a forced flow of dry air. Thus the efficiency of the drying peocedure affected the MIT as is also the case with milk powders. Thus damp powders are much more likely to ignite than correctly dried powder (Synnott et al., 1986). Conventional reaction kinetics can not describe the reaction between dry pie wool and oxygen as it includes a complex reaction involving the oxidation of an olefin. Walker et al. (1982) uses a single temperature coefficient for the reaction rate and reaction time in constructing a plot of -55- equivalent reaction rate at 0°C versus equivalent time at O‘C. Walker et al. (1965B) examines the particular problems of applying ignition theory to porous solids. Knowledge of the kinetic data is a prequisite for such a study. Due to varying temperatures across a reacting sample, reaction rates also varie. Using a calculation developed previously ( Walker et al., 1965A), the reaction rate values were determined. using' 'mean' temperature ‘values, taken from experimental time / centerpoint temperature profiles of cylindrical piles of acetone extracted dry wool. The 'k' values are found to decrease ’with time at constant temperatures. The authors conclude that the normal understanding of order of reaction and of Activation Energy does not apply to the reaction between dry wool and oxygen. Problems associated with a)the heterogenous nature of the gas/solid reaction, and b)the speed of the main (explosive) reaction, place the reaction in a different category from the conventional solid state reaction and associated theory. An equation is reported by Walker et al.(1969) for calculation of the ignition temperature of porous solids, with the surface temperature close to ambient. While this restricts the general application of the model, it does have the advantage of not directly involving reaction rate -55- in the equation. The effects of sample size and bulk density on the ignition temperature of a cylindrical pile of scoured wool are determined experimentally . The authors dismiss as unlikely the zero order nature of the exothermic reaction of solids with air, in spite of some results obtained by other investigators (e.g. Thomas et al. , 1961A) . WM Due to the growing interest in curtailing dust fire / explosion phenomena in the dairy industry, a catalogue of milk powder ignition data has been built up in recent years. Duane et al. (1981) published the results of a study on various milk powders including whole and skim milk, whey powders, and powders containing 22-30% tallow or 26% coconut. The powder sample were suspended in wire mesh cubes and heated until ignition ( or non-ignition) occurred. The Minimum Ignition Temperatures (MITs) were thus determined for the samples (see Table 3.1). Sample size, packing density and particle size were varied to assess the effects on the MIT values. The authors conclude that fat-filled powders are not inherently more likely to ignite than conventional powders. -57— 13h11_1&1 MIT data for milk powder cubes (Duane et al.,1981). man an: 30% Tallow 187'C 22% Tallow 188°C Skim milk 182'C Buttermilk powder tested soon after manufacture had an MIT of 170'C, whereas a similar sample, tested after protracted storage had a higher MIT value of 177'C. This conclusion is important in pinpointing critical stages in the powder production cycle. Beever (1984) used the standard technique as outlined in section 3.1.2.2 to determine values of MIT for cubes of milk powder with sides of 50, 75 and 100 mm. Taking a value of 2.64 for 5cr for a cube and 0.92 and 3.46 for a layer and sphere, respectively, the author used the results to determine the critical layer and sphere sizes for the powders tested. Ignition in dust layers is clearly a problem in spray dryers as there is a tendency for layers of dust/powder to build up on the dryer walls or in the crevices or corners. The regions where this problem may occur are indicated in Fig. 3.3,. This is more likely to happen with damp or 'sticky' powder or on a damp surface. When such clumps form, they often break loose and give rise to the so- called 'snowball' effect (Synnott et al., 1986). The ability to predict a critical sphere diameter is thus of value in determining the likelihood of this effect giving rise to an ignition situation. The results in Table 3.2 -53- Hot Air Inlet Duct Atomizer Fine: Return Pipe Exhaust Air Duct Fig 3.3 : Typical Conical Sprag Drier Section Showing Temperature Ranges and Powder deposit Areas. -59- were obtained by Beever (1984) using cubes with sides of 75 mm. 1§h1g_1;z MIT values for milk powder cubes (Beever,1984). REEDEB HIT Skim milk 156.0'C Skim + 22% Coconut oil 147.5'C Skim + 30% Tallow 149.5'C Skim + 30% Tallow + 10% Soya 148.0°C The difference between the skim and the fat-filled powders is about the same as that noted by Duane et al. (1981). The MIT values given by Beever(1984) tend to be lower than the corresponding values of Duane (1981). O'Mahony et al. (1988) determined the MITs of a 30% tallow fat-filled milk powder for various shapes including cubes, spheres, _ cylinders and layers. They investigated how the standard equation, Equation [ 3.33 ], is used to predict criticality results for the main regular geometries based on tests conducted with one shape. Their conclusion was that extrapolation from spherical test results gives the best correlation with experimental tests performed on layers on a hotplate. Estimating critical dimensions such as safe silo sizes is thus best done using spherical test samples. The simulation model developed in the present study also considers spheres. For the powder tested by O'Mahony et al. MIT values varied from 131'C for a 51 mm -50- sphere (radius) to 206'C for a 5 mm layer of powder. In general, increasing the sample size gives a decreasing value of the MIT as the ability to dissipate the heat decreases as the thermal resistance increases. This is what is expected intuitively. It is also predicted by the approximate model of Boddington et al. (1981). Synnott et al. (1986) studied fat-filled powders to quantify the effect of fat percentage on the MIT. Results show that the MIT tends to decrease with an increase in fat percentage. Unsaturated fats are more prone to oxidation than saturated fats. By checking the level of unsaturation in a number of commercially manufactured fat- filled powders - 26% coconut, 26% palm kernel, 33% lard and 44% butterfat - it was found that the MIT decreases with an increasing total of unsaturation. Differential Thermal Analysis curves were obtained for a number of fat- filled powders to identify the main 'thermal' reactions involved in heating/combustion. Using a cross-flow oven, it was found that the MIT decreases with an increas in the air velocity. Prompted by a minor industrial fire, caustic soda was added to a powder formula to simulate contamination of dryer feedstock with CIP (Cleaning in Place) detergent. The MIT was found to be reduced considerably (e.g. from 176.5'C to 160.5'C in one sample). Amounts of less than 1% added caustic were sufficient to decrease the MIT significantly. Rewetting of powder, -51- another common occurrence industrially, was also found to give a dangerous reduction in the MIT value. 11111_IDI_§!HEABI_112§11 The International Dairy Federation Expert Group on spray drying of milk powder summarised the self-ignition risk as follows (IDF, 1987): 1) All dried dairy powders can proceed to self- ignition. 2) Powder sample dimension has a substantial influence on ignition: the MIT decreases by 15K as the sample radius doubles. 3) The milk powder fat content does not significantly affect combustibility. 4) The smaller the powder particle size, the greater the fire risk. 5) The low-temperature crystallisation of lactose can sometimes lower the MIT value of the powder. 6) With a sample critical dimension up to 50 mm the MIT is above 130‘C: larger clumps or layers ( up to 150 mm) can self-ignite at 80- 90'C. 3I5_SALQBIHIIBI_AHD_SIEIIIQ_EIEDI£§ One of the key factors in determining whether an oxidation -52- reaction of the type under study proceeds to combustion or explosion is the rate of heat generation. Equation [3.2], repeated. below, summarises the key terms in the heat generation phenomenon : c - Q cin f e'EA/RT {3.2} The value of each of the variables plays a significant role in determining whether or not the exothermic oxidation reaction will proceed to a hazardous ignition or explosion. To accurately model the reaction, each of the terms must either be determined experimentally or otherwise derived from the ‘literature. Typically, simplifying assumptions are made about the main energy generating reaction to allow the mathematical analysis to proceed. Alternatively, composite values of variables are employed. 31111.8EAQIIQE_BAIE§ Walker et al. (1977) conducted experiments, involving the controlled release of electric heat, to evaluate a calorimeter equation which allows the calculation of the reaction rate from time / temperature profiles of the centerpoint of a sphere of wool with a perfect thermal contact with the ambient( i.e. negligible surface thermal resistance) .' The thermal conductivity and specific heat -53- were determined experimentally and modified to take account of the experimental conditions for the ignition tests. When the wool was heated at a constant rate, the heat generation rate was evaluated to within 7%. Tested at diminishing reaction rates, however, the accuracy decreased as the time index increased. Greater time dependence led to larger error in the reaction rate calculation, increasing from about 10% for a relatively low dependence to 30% when time played a more important role in the rate of reaction. An eigenfunction transform method was proposed by walker et al.(1978C) to handle a reaction rate which diminishes as a power function of time, for a self-heating sphere with. the Dirichlet boundary’ condition. This allowed a reaction rate to be calculated from an observed central temperature rise. For reaction rates up to 0.3 the error in the heat output calculations was of the order of 10%. At a rate of 0.7 the error increased to 40%. The model did offer an improvement when dealing with two simultaneous but independent reactions occurring in the test material. Walker's preoccupation with centerpoint profiles and calculations loses some significance with the milk powder samples under study as the sphere center may not always be the center of ignition. Walker et al: (19788) used a finite difference technique -64.- to calculate the central temperature rise in a reacting sphere where the rate increases as an exponential function of temperature and decreases as a power function of time. The calorimeter equation. of ‘Walker' et al. (1977) was evaluated under the influence of a temperature coefficient, with the surface temperature maintained at ambient. Heat generation was written in terms of an equation which does not include EA, the Activation Energy. Thus, using standard notation, z = Ka[2]KbT [t]Kc [3.42] where Z - Reaction rate; Ka, Kb and Kc - Reaction rate constants: t - temperature: T - time (s). Allowance is made in the model calculations for the fact that the exothermic reaction has progressed a finite amount prior to commencement of the simulation calculations. They note the statement of Hinshelwood in 1929 that a marked deviation from linearity in a semilog plot. of reaction. rate [against ‘the reciprocal of temperature is evidence that the reaction under investigation is made up of at least two concurrent reactions. The rate of oxygen uptake is possibly a better indicator of the reaction rate than the rate of heat generation. Walker et al. (1983) use the equation again to -55- calculate the central temperature rises in reacting spheres (with the surfaces maintained at ambient). The conventional rate laws are based on a fraction of reactant remaining and include an Activation Energy term. These laws do not account for the oxidation by two (or more) simultaneous reactions. A. simple temperature coefficient of reaction is a more useful and accurate concept. Inaccuracies in assigning values to EA and the pre-exponential factor in the Arrhenius equation give rise to an error-compensating effect. Varying degrees of success and failure are quoted where workers use Equation [3.42] to model self-heating reactions. The authors show how both temperature coefficient and rate constants can be derived from calculated centerpoint temperature profiles. The same profiles can be used to assess the effect of time and temperature on the reaction rate of reactions which fall off with time and increase in rate with increased temperature. An approximate analytical model has been proposed by Boddington et al. (1980) to describe an exothermic reaction in a reacting sphere with the surface at ambient temperature and with a diminishing reaction rate. Finite differences were used together with a time-varying value of the 6 parameter to predict time / temperature profiles, times to ignition and other criticality data: the standard -66- Frank-Kamenetskii assumptions were made along with the time dependence of the reaction rate. The results include a useful family of curves illustrating the various possibilities of the critical and non-critical conditions for the F-K heat transfer model. The Differential Scanning Calorimeter (DSC) is used to monitor energy changes in systems undergoing a thermal process ( Moshenin, 1980 ). Its widest application is in determining specific heat and it has been employed in cataloguing many food products in this regard. It has also found application in trying to determine other thermal properties of foodstuffs. Thus, Lovric et al. (1987) employ the DSC to determine factors such as heats of melting and fusion as well as specific heat for a number of liquid and semi-liquid foodstuffs. Quinn et al. (1980) use the DSC to assess changes in heat stability of meat protein during processing of meat into sausage batter. A Calvet type of heat flow calorimeter was used by Raemy et al. (1982A) to determine the specific heat of coffee and chicory in order to understand the exothermic, self-heating reactions of these products. The heat absorotion was measured simultaneously in an empty sample cell and in a cell with product, as both samples underwent the same temperature rise. In general, the specific heat -67- increased as temperature increased. The results reported are mainly for dry products: the authors suggest that small moisture differences have a considerable influence on the measured specific heat ( Cp ) values. A heating rate of l'C/min was used to the detect exothermic reactions. ( due to roasting’ and carbonisation of the samples). Further' tests also identified. exothermic peaks in the heating of wheat and whole rice. Raemy et al. (19828) catalogue the Reaction Enthalpies and threshold temperatures for a number of cereals such as wheat, maize and rice using Differential Thermal Analysis. As the samples were heated in sealed pans at a controlled rate (1 ‘C/min ), from ambient up to 270'C, they underwent exothermic reactions. The enthaply of these reactions are shown to be closely related to the carbohydrate content. Raemy et al.(1983) outline the use of Differential Thermal Analysis (DTA) and a Calvet type Heat Flow Calorimeter in determining the exothermic behaviour of selected carbohydrates. While the Calvet technique does not permit detailed analysis of the individual decomposition reactions, it does provide a useful record of the overall process taking place when carbohydrates are heated. The behaviour of food components is studied by Raemy et al.(1985A) using both DSC and DTA. To simulate industrial conditions such as freezing and roasting , for example, the instruments had to be operated outside their normal operating range. In an extension to this work, Raemy et al.(19858) reports on the use of high pressure DTA tests to assess self-ignition properties of food powders. Raemy et al.(1985C) adapts thermal analysis techniques to study the dust explosion phenomenon W The standard generalised rate of reaction is (Widmann, 1982) do /dt = k(l-c)n [3.43] where dc/dt - rate of reaction, 3'1 k - reaction rate constant, 5"1 a = degree of conversion or fraction reacted ( a - 0 at t - 0 ) n = order of reaction The Arrhenius equation relates the reaction rate constant to temperature as : k = k0 exp ( -EA/RT ) [3.44] where k0 is the pre-exponential factor, and EA, R and T -59- are as previously defined. Thus, da/dt a k0 e'EA/RT (1-a)n [3.45] Each incremental reactant component, do, produces a corresponding enthalpy change dH (i.e. the fractional incremental degree of conversion is the same as the fractional change in enthalpy ) : da = dB / AHTOTAL [3-46] Taking a time derivative : sis-iii - _1_ [3.47] dt dt AHTOTAL Thus, the rate of reaction is directly proportional to the DSC signal (i.e. power input). Equation [3.46] may be written as: a = mPART / AH'ro'I'AL [3-48] Thus, the degree of conversion is proportional to the associated enthalpy change. If AHr is the remainder of the enthalpy curve, (Figure 3.4), -70- AH AH, Fig. 3.4 Dsc Enthalpy Curve. -71- 1 " a = AHr / AHTOTAL [3-49] Rewriting [3.45] in terms of enthalpy yields : on . 1 = k0 e‘E-A/RT . AHr n [3.50] dt MTOTAL AHTo'rAL Taking natural logarithms gives : ln H' = ln k0 - EA 4. n In AH:- [3.51] AHTOTAL RT AH'ro'I'AL where H' = dH/dt. The DSC records H' and T throughout the experiment, and performs the integrations AHr and AHTOTAL as instructed by the preset program at the conclusion of the run. Three unknowns are calculated from the H vs T curve ( i.e. k0, EA and n ). A multiple linear-regression analysis is performed to calculate the unknowns. They are printed out with confidence limits computed on the basis of 95% probability (see Figure 4.6). -72- MW Differential Scanning Calorimetry is a technique whereby the physical and chemical properties of a sample may be examined when subjected to a defined temperature program. The temperature program ensures that the sample is measured at a constant temperature (isothermal program) or with a linearly-increasing temperature (dynamic program). A graph of the enthalpy change vs. temperature/time may then be examined. Physical transitions may include fusion, re- crystallisation, evaporation, sublimation, condensation, solid-solid transition, glass transition. Chemical transitions include thermally induced decomposition, oxidative decomposition, polymerisation, polycondensation and specific heat (Widmann, 1982). Thus, significant information about a sample can be determined by Differential Scanning Calorimetry. Using Differential Scanning Calorimetry, the temperature of the sample is compared with the temperature of an inert sample (air in this case). The temperature changes which occur during the physical or chemical changes are detected by a differential method. The advantage of the DSC technique over thermal analysis is that the temperature of the sample, Ts, is recorded in thermal analysis as a -73- function of time, and. a heating or cooling curve is recorded: small temperature changes occurring in the sample are generally detected by this method. With DSC the detection thermocouples measure differences between the sample temperature (Ts) and the reference temperature (Tr) . Where these differences are small they can be detected with an appropriate voltage amplification device. It also allows the use of very small samples (mgs). Graphs of the two techniques are shown in Figure 3.5. The essential difference between the two curves in Figure 3.5 is that in TA no enthalpic transition is monitored, while in the DSC analysis exothermic and endothermic changes occur. Since no other temperature changes take place in the sample undergoing Thermal Analysis, no deviation from the linear temperature is detected in the ‘ sample temperature. However, in the DSC deviations occur at the programmed initial reaction temperature, T1, due to the temperature changes caused by the exothermic or endothermic reaction. These changes are computed with respect to Tf, the final temperature,and the temperature of the sample returns to that of the system. From the DSC graph the difference in temperature (Ts - Tr) is recorded as a function of the system temperature, T. At T1, the curve deviates from the horizontal position to form a maximum or minimum peak, depending on the enthalpic change. -74- DSC D T A 4 u. 4 p- .- Exothennic t-‘o' "‘ - ' Ti V Endothenmc a [- e '1‘: T1 1} Templ‘l'ime Tempmme Fig. 3.5 DTA and DSC Compared. -75- The completion of the reaction temperature, Tf, does not occur at the curve maximum or minimum but rather at the high-temperature side of the peak. The exact position depends upon the instrument arrangement. Thus, in the differential method, small temperature changes can be detected : also the peak area is proportional to the enthalpic change and the sample mass. The size, shape and position of the peaks yield different information about the sample, and can be used for qualitative identification of the sample material. Also, as the area under the curve is jproportional to the heat change, the 'technique is useful for the semi-quantitative or, in some cases, quantitive determination of the heat of reaction. Thus, as the heat of reaction is proportional to the amount of reacting substance, DSC can be used to evaluate quantitavely the amount of substance present if the heat of reaction is known. Hence, the technique finds extensive use in the qualitative and semiquantitive identification of organic and inorganic compounds such as metals, minerals, fats, oils etc. Quantitatively, it can be used for the determination of a reactive component in a mixture, or the heat of reaction in physical and chemical changes (Raemy et al., 1983). ii§i§_§£l§IIIQ_EIAI To determine specific heat values for the milk powder -75- samples, the predictive equation developed by Heldman et al. (1981) is used in conjunction with the detailed compositional data. This equation allows a composite specific heat value, Cp, to be determined from : cp = 1.424 xc + 1.549 xp + 1.675 xF + 0.837 xA + 4.187 xx [ 3.52 ] where X refers to the mass fraction of the various components of the material i.e. C : Carbohydrate: P : Protein: F: Fat: A : Ash: M : Moisture. The numerical coefficients are the respective specific heats of these components, 1.675 kJ/kg K being the specific heat of the solid fat phase. 21111.22EEIIX Milk powder density plays an important role in this analysis. This is calculated using the compositional data by' adding ‘together 'the [component densities as follows (Heldman et al., 1981): .l-l‘n+_xr +1211: P PM PF Psur [ 3.52 ] The subscripts are as previously defined. SNF is the Solids Non Fat component. Carr (1976) quotes a value of 610 kg/m3 for milk powder density. -77- The Finite Element method is generally associated with topics of civil engineering analysis and design, since its original application was in the area of structures. Davies (1980) traces the first attempts to develop the analogy between discrete elements (such as bars and beams) and the corresponding sections of continuous solids to Hrenikoff (1941) and McHenry (1943). Clough (1960) was the first to introduce the actual term 'finite element' to this form of analysis. In the mid-60's the Finite Element analysis was extended to dynamic problems. Extension to non-structural studies followed including the transient heat conduction problems (Wilson et al., 1966). As outlined in more detail below, Anderson et al. (1974) and Misra et al. (1979) have also looked at problems of heat transfer. Segerlind (1984) gives a comprehensive treatment of the general Finite Element application, including the problem of convective heat transfer. -73- C L SIMULA ION OP SE - N IQEIIIQH Wake et al. (1973) outlined a variational technique to solve a non-linear eigenvalue problem, (e.g. the steady- state thermal ignition problem). The method allows calculation of a theoretical critical parameter analogous to the Frank-Kamenetskii critical parameter, 5cr- A theoretical study for a first order reaction occurring due to assymetric heating of a slab was published by Tyler et al. (1981). They used a finite difference technique to predict the temperature and reactant concentration as functions of time and location. The exposed slab face was subjected to Newtonian cooling and the slab's properties,including the thermal conductivity and specific heat, were assumed independent of temperature. Kordylewski (1980) studied a complex reaction for a porous body in when the heat generation is first order with respect to both the porous solid fuel and oxygen concentration. He concluded that the critical parameter, 5cro depends on the ratio of the Lewis number (Le) to the dimensionless adiabatic temperature rise. The Lewis number is defined as the ratio of the Prandtl number for heat transfer to the Schmidt number for diffusion : Le Pr / Sc , [3.54] -79- where Pr = kinematic viscosity / thermal diffusivity and Sc 8 kinematic viscosity / diffusion coefficient. Shouman et al.(1975) studied the onset of thermal ignition in a reactive slab with unsymmetric boundary temperatures (i.e. one side heating, the other side cooling) under steady state state conditions. They concluded that the standard 6 parameter was of little use in analysing the unsymmetric case. Boddington et al. (1981) developed an approximate, analytical solution for spherical reactants in. a non— isothermal steady state using the standard Dirichlet boundary condition. To account for the effects of self- heating on the reaction rate parameters, they defined a correction factor, v: v = Q R A exp(-EA/R Tamb)/(K+l)k(dT/dr)r=ro [3.55] where the variables follow standard notation as defined previously. The presence of the exponential term means a numerical solution must be found. A quintic approximation was proposed whereby exp(-EA/RTamb) was replaced by a fifth order polynomial in e, the F-K dimensionless temperature excess (cf. Equation [3.25]). Alternatively, a second order reversion was used to approximate an infinite series solution to the heat balance equation at criticality. The resultant surface temperature profile was -80- used to evaluate v. For the case of an arbitrary Biot number, the approximating equations were recast to express 8i in terms of 6, the F-K parameter. Though the model is mainly structured in terms of the traditional central temperature excess, the prediction of a varying maximum temperature position is also possible. Accuracy to within 2% is claimed when results are compared with known exact solutions. The authors proposed an extension of the quintic approximation to the temperature dependence of thermal conductivity, depicted as an exponential by Wake (1980) . A finite element approach .to the problem of spontaneous ignition in a theoretical reactive solid has been presented by Anderson et al. (1974) . Their study centers around determining values of the Frank-Kamenetskii critical parameter, 6cr, for general shapes of samples undergoing a zero-order exothermic reaction. Boundary conditions are of the Dirichlet type, with one example of a finite surface heat transfer coefficient also covered (the case of a cylinder in steady state) . The discretised finite element equations were established using the Galerkin's method which were solved by an incremental procedure. Sample properties were assumed invariant with respect to temperature changes. In solving the transient problem, theoretical values of material properties were used to demonstrate the scope of the model for determining -31- 5cr for a sphere. Assumed 'lumped' values of kinetic data were used throughout. Details of the shape functions employed and some of the model parameter values are not given in the publication. Working with an ambient temperature of 500K, thermal and reaction parameters were chosen to give values of 6 of 2,4,8,16 and 40, respectively. The resultant induction times were calculated for a sphere. 5cr was found to be 3.32, the value below which ignition does not occur. Taking different values of 6, sphere temperature profiles were plotted at various ambient temperatures. The profiles show that as 6 increases the nucleus of ignition moves away from the sphere center. The sphere modelling work assumes F-K boundary conditions. No experimental verification of the simulation was attempted here. -82- 1i2_IIRIBIEEEIAL_IEIBQDEQIIQE The primary tool developed in the course of this work is the mathematical model to be used for sensitivity studies on self-heating and self-ignition of commercial milk powders. A laboratory phase was essential to the project for two reasons. Firstly, the model, as it is designed, needs seed data from classical oven tests or from DSC studies to initiate the simulation process. Information on the experimental conditions is also needed for the model. Secondly, the time / temperature profiles produced by the model need experimental verification. While an elementary version of the Finite Element model can be verified (using an analytic solution for the simple case of heat transfer with. an inert sample :material), no analytic benchmark exists for the case of a product undergoing an exothermic reaction. Thus, an experimental verification is necessary for the model. The oven test to decide Minimum Ignition Temperatures for different samples / sample sizes often involves lengthy test runs, particularly to identify the 'lowest non- ignition' temperature ( sometimes over 2000 mins ). The Differential Scanning Calorimeter offers a faster and potentially more accurate technique to obtain the essential kinetic parameters necessary to run the -83- modelling program. In the experimental phase of the work milk powder samples are analysed using the DSC. ii1_IIIIBIHIEIAL_2ABAHITIB_BAEQI The range of experimental parameters used in this phase of the study was decided by two principal factors : 1)The conditions likely to be encountered in industrial milk-powder drying or storage facilities. 2)The capability of the experimental apparatus to recreate and monitor accurately the industrial conditions. As with all laboratory-based studies, scale-up factors pose problems. In this study one such problem manifests itself clearly. This is the type of limitation imposed on conclusions drawn from experimental results due to the necessarily limited size of a powder sample which can be studied in a test oven. Reference has already been made to the significance of the reactant consumption in understanding self-heating and exothermic reactions ( cf. Section 3.2.2 ). As the reactant is consumed, less reacting matter is available to continue the 'runaway reaction' .Thus, the reaction may not in fact reach an ignition situation. This may also happen industrially when -84- the subcritical self-heating occurs in an isolated clump of powder, resulting in. a charred. mass in the dryer chamber with no attendant fire or explosion. In this case, oven tests are very close to the real situation. Experimental results appear less effective, however, when trying to extrapolate results to industrial situations where, once a reaction enters the self-heating spiral, it will have a ready supply of reactant (i.e. powder and oxygen) to sustain the reaction until critical ignition occurs e The traditional theory of Frank-Kamenetskii views the system as one of unlimited supply of reactant. An advantage of the simulation in this study over the experimental work is that it may be extended, rather than extrapolated, to mirror large-scale industrial conditions. While not immediately accounting for reactant consumption, this latter condition becomes prominent in the post- ignition situation when reactant consumption increases dramatically with exponentially increasing temperature. In the present study, where the intention is to avoid criticality. Winn ”WWW W As the present work is concerned with simulating industrial conditions, the experimental temperature range -35- is limited to the typical temperature range found in conventional spray-drying situations. The range of temperatures is shown in Figure 3.3 for a typical conical spray dryer (Duane et al. 1981) . The temperature ranges from 40°C at the product inlet, the atomiser, to 200°C at the hot air inlet. The dryer surface temperature varies from 80°C up to 100°C. If a further source of heat develops (e.g. external surface welding, radiation from a spot lamp, etc.) , the surface temperature may increase significantly. As the diagram also shows, product deposits may accumulate in angles or corners in the dryer chamber, inhibiting heat dissipation and giving rise to increased critical dimension, and, as will be shown later, increased risk of self-ignition. Thus, the oven temperature range investigated is between 127°C and 380°C, encompassing the span of temperatures encountered in a commercial dryer. W The oven employed to conduct the heating / ignition experiments is a modified Townson & Mercer convection oven. Additional heating capacity was built into the oven by including extra heating elements. The fan speed in the oven is controlled using a 0-240 V, 2A variable alternating-current transformer to power the fan. This allows variation of air speed and hence the surface heat transfer coefficient. The original oven thermostat was -35- replaced ( to provide accurate control of the oven temperature )with a Honeywell CL-40 three-term controller capable of maintaining the temperature within 11°C over a temperature range of 0 - :300'C: it employs a .J-type Copper/Constantan thermocouple. The controller includes a power compensation circuit which maintains a constant power output in case the line voltage changes from -15% to +10% (which. are the typical supply 'variation extremes within the laboratory). Steady temperature is achieved rapidly (within 10 minutes from start-up) and reliably using a controller. Figure 4.1 shows the basic layout of the test oven. The test sphere is suspended in the oven, as shown, with thermocouples used to monitor and record the temperature within the sphere and in the oven itself. The J-type thermocouples give a measurement accuracy of better than 11°C. The thermocouple voltage is measured on a Prema 5000 DMM/Scanner digital multimeter: the time / temperature information is communicated to an Apple computer via an IEEE-488 interface card. Alternativery, a Philips multi- point recorder is used to record the system temperature. A shrouded thermocouple is employed to measure the ambient temperature in the oven. This is necessary to minimise possible errors due to radiation. -87- Oven Temperature 1— Councilor Multi- Channel Temperature . . . Recorder Spherical Sample Lmlbe I ‘ 4‘ Fan Speed Control Fig. 4.1 Test Oven. —88- .!!!"x 1 O O ”"191-51'1. ,!i__1_=."; 0?" iii! A relatively simple but effective experiment was used to determine the effective surface heat transfer coefficient for given conditions in the oven. Figure 4.2 shows a typical time / temperature plot for an aluminium sphere. The sphere was suspended in the air stream and, using thermocouples, the centerpoint temperature was recorded. As aluminium is an excellent conductor of heat, the major source of thermal resistance within the sphere is located at the surface due to the prevailing boundary condition. Quantification of the boundary condition is required for the simulation. An equation for heat transfer in a sphere with negligible internal but finite surface resistance to heat transfer is used (Wong, 1977): Tinit ’ Tamb [4'1] where Tinit is the initial temperature of the sphere and the other variables are as previously defined. All parameter values are known except the heat transfer coefficient, h, and the temperature / time data , T and t, respectively. The value of 'h' is calculated at points taken at regular intervals along the curve, allowing for the initial, short time lag across the sphere radius. The -89- 120 r 100 '- 80*- ('C) Temp. 60'- 40" l a l 20 a J A l A l A I a 0 10 20 30 40 50 60 Time (mins) Fig. 4.2 4" Aluminium Sphere Heat-up Curve. -90- average 'h' value is calculated from these values. For a given fan-speed and sphere size, 'h' values do not vary greatly over the span of the heat-up curve. Details of the calculations are shown in Section 6.1. 1121211.!ABIIQL§_D§E§IIX A Beckman Model 930 Air Comparison Pycnometer was used to establish the particle density of the powder samples. A sketch of this is shown in Fig. 4.3. The pycnometer consists of two chambers, a reference chamber and a sample chamber. The powder is initially weighed and the difference in volume due to the powder is noted. A simple calculation then gives the particle density. mm A jolting volumeter was employed to measure the apparent or bulk density of the powder samples. This apparatus consists of a platform mounted on a shaft supported in a vertical position by a sleeve around the shaft, which in turn rests on a cam wheel (see Figure 4.4) . The powder sample is placed on the platform in a graduated cylinder. -91- 4* Dmtfferenm Reference meme: Piston Measuring Piston. fi 23!!) I Strung = : Tara Number Sample 351— d—-)§ V x CUP (in cc) Fig. 4.3 Pycnometer. -92- Graduald Stand cylinder Gearbox Counter Moor Conner rem UDDUDU Poveronloff \ ...-I- g L ‘1 I fi— Fig. 4.4 : Jolting Volumeter. -93— The cam is rotated by a constant-speed motor at 250 rpm. On each revolution the platform is gradually raised and dropped through a distance of 3 +/- 0.1 mm: 1250 revolutions are allowed before the volume of the sample is measured. Knowing the weight of the powder, the apparent density can be calculated ( Foley et al., 1974). The dairy-powder samples used in the experimental phase were commercially-produced powders (i.e. manufactured under industrial conditions in a commercial milk powder plant). To [ensure ‘uniform lheat transfer conditions across all samples the sample bulk density is carefully controlled for the ignition tests. The powder bulk density is -94- measured using the jolting volumeter as previously outlined. Steel mesh spheres are constructed to hold the powder in the convection oven. Knowing the exact radius of the sphere, the correct weight of powder is determined to ensure the sample is at the correct bulk density prior to beginning the heat transfer experiment. This calculated weight of powder is measured out and put into the sphere, tapping the sphere until the full weight of powder is enclosed in the sphere. The samples were subjected to a convective air-stream in the stainless steel mesh spheres. The spheres were constructed from a No. 30 gauge mesh to the required dimensions. As the steel mesh is of high thermal conductivity, it is assumed. to present. no significant thermal resistance to the transfer of heat from the powder surface to the ambient air stream. Actual dimensions of the spheres are recorded for use in the calculations. The nominally 4" sphere has an actual measured diameter of 3.974", the 3" sphere an actual diameter of 3.166", etc. Table 4.1 lists the nominal and actual diameter and also shows the calculated 'h' values. -95- xgplg_1;1 Powder sample radii and associated 'h' values. flgminal Diameter Actual Diameter e 've ' 'v (inches) (inches) W/mzK 4 3.974 19.21 3 3.166 20.12 2 2.148 21.26 1.5 1.688 21.78 ii1_IIH3.1.123223AIEBI_EBQIIL£§ To determine the combustion parameters, the traditional oven procedure was employed. Using a thermocouple located at the sphere center, the time / temperature profile was recorded for' a number of different spheres. For each sphere the oven temperature was varied to find the lowest temperature at which a sample of particular radius ignited, and the highest temperature at which the sample failed to ignite. The initial temperature used in the oven is chosen on the high side of the assumed MIT. Known kinetic data, derived either from DSC tests or from oven tests on similar powder samples, may be used in the model to determine approximate MIT values around which to base the experimental design. This makes a significant reduction in actual experimental work. Typically three or four oven runs are required before a non-ignition temperature is encountered. For example, Figure 4.5 shows three curves for a 4" sphere: 138.9'C is the lowest -96- ('C) Temp. mm P a l 0 1000 2000 Time (mins) Fig. 4.5 4" Milk Powder Ignition Curves. -97- ignition ‘temperature and 137.27'C is the .highest non- ignition ‘temperature. The time ./ 'temperature. data for these curves are listed in Table 4.2. In the case of the highest non-ignition temperature, the experiment was allowed to continue for 2000 mins ( :> 33 hours) before. discontinuing' the recording. In. the IDF Bulletin (IDF, 1987), 24 hours (i.e. 1440 mins) is the upper time limit recommended. In the case of the other two ignition curves, the sample is seen to proceed to a runaway exothermic reaction (self-ignition), reaching a: temperature peak, after which the sample temperature decreases due to the finite limitation of reactant material. The recording ceases at this point. The MIT for the powder in Figure 4.5 is 138.08'C, obtained by taking an average of 138.9°C and 137.27'C. A 8 N6 CA R ER A thermo-analytical technique was employed to identify the kinetic data characterising the exothermic reaction in the self-heating milk-powder samples. A Differential Scanning Calorimeter (DSC) was used in this phase of the study. The particular model was the Mettler DSC 20 . _98- 13211—431 Time/Temperature data for 4" sphere. Amb.137°c Amb.139°c Amb.l43°c Time 19221 Tenn; Tease (mins) ('C) ('C) ('C) 2 20.15 20.51 19.64 6 17.85 18.30 17.91 10 18.15 18.50 18.47 14 19.75 19.62 20.22 18 22.82 22.09 23.38 22 26.94 25.79 27.61 26 31.49 30.16 32.42 30 36.00 34.71 37.22 34 40.21 41.02 43.82 38 44.02 42.94 45.88 42 47.57 46.58 49.80 46 51.09 50.08 53.78 50 54.74 53.73 57.88 54 58.28 57.53 61.84 58 61.59 61.13 65.46 62 64.59 64.40 68.63 66 67.26 67.28 71.34 70 69.57 69.75 73.63 74 71.57 71.85 75.56 78 73.30 73.65 77.19 82 74.79 75.17 78.59 86 76.07 76.49 79.80 90 77.20 77.62 80.84 94 78.20 78.59 81.78 98 79.09 79.44 82.63 102 79.89 80.20 83.40 106 80.62 80.88 84.14 110 81.31 81.51 85.16 114 81.95 82.10 85.81 118 82.57 82.66 86.44 122 83.19 83.20 87.02 126 83.80 83.74 87.59 130 84.41 84.29 88.17 134 85.03 84.84 88.77 138 85.65 85.40 89.43 142 86.29 86.00 90.17 146 86.95 86.61 90.99 150 87.64 87.29 91.88 154 88.34 88.00 92.86 158 89.07 88.77 94.47 162 89.83 90.03 95.62 166 90.62 90.92 96.85 170 91.43 91.89 98.11 —99- 92.29 93.18 94.08 95.04 96.00 96.98 97.99 99.00 100.03 101.07 102.09 103.11 104.13 105.14 106.12 107.10 108.05 108.99 109.90 110.80 111.67 112.52 113.38 114.20 115.00 115.77 116.53 117.26 117.98 118.67 119.33 119.97 120.58 121.19 121.77 122.33 122.90 123.44 123.99 124.54 125.07 125.60 126.14 126.68 Alb.139'C TOME, ('C) 92.90 93.97 95.07 96.21 97.38 98.57 99.77 100.99 102.19 103.41 104.61 105.79 106.98 108.14 109.27 110.39 111.48 112.55 113.58 114.59 115.57 116.52 117.44 118.34 119.21 120.06 120.88 121.68 122.48 123.27 124.45 125.23 126.02 126.83 127.65 128.48 129.32 130.16 130.99 131.41 132.23 133.02 133.78 134.52 -100- Amb.143°C 1:22; ('C) 99.42 100.77 102.14 103.52 104.93 106.34 107.75 109.17 110.57 111.93 113.27 114.58 115.86 117.10 118.31 119.50 120.66 121.80 122.93 124.04 125.16 126.28 127.42 128.58 129.75 130.93 132.13 133.29 134.39 135.47 136.50 137.51 138.52 139.54 140.60 141.68 142.81 143.98 145.21 145.85 147.16 148.54 149.99 151.52 Tine (mins) 350 354 358 362 366 370 374 378 382 386 390 394 398 402 406 410 414 418 422 426 430 434 438 442 446 450 454 458 462 466 470 474 478 482 486 490 494 498 502 506 510 514 518 522 Tabls_ilz_lsentldl Amb.137°C Temp, ('C) 127.22 127.76 128.30 128.82 129.33 129.85 130.33 130.80 131.27 131.72 132.15 132.58 133.00 133.39 133.79 134.19 134.57 134.97 135.36 135.76 136.16 136.55 136.95 137.35 137.74 138.14 138.54 138.94 139.33 139.74 140.15 140.97 141.38 141.80 142.22 142.65 143.07 143.50 143.94 144.38 144.83 145.29 145.52 146.10 Alb.139'C ISERe ('C) 135.24 135.59 136.30 137.00 137.73 138.45 139.19 139.19 139.95 140.73 141.51 142.32 143.15 144.01 144.90 145.83 146.78 147.29 148.31 149.93 151.09 152.31 153.62 154.99 156.49 158.12 159.93 161.98 164.34 167.05 170.06 173.29 176.67 180.13 183.62 187.15 190.81 194.80 198.80 203.32 208.30 213.78 219.80 226.27 -101- Amb.143°C ISEEs ('C) 153.12 153.96 155.71 157.56 159.57 161.78 164.32 164.32 167.32 170.91 175.06 179.67 184.55 189.51 194.59 199.93 205.62 211.73 218.27 225.16 232.28 239.51 246.72 253.70 259.99 265.20 269.37 272.65 275.34 278.73 280.88 283.17 285.91 289.39 293.71 298.82 304.69 311.12 318.01 325.49 333.76 342.78 352.58 363.32 In the DSC, a sample is heated at a controlled rate. The heat input is controlled to allow the sample temperature to increase linearly. By monitoring the heating rate and temperature, the progress of the reaction is recorded as shown in Fig. 4.6. This printout has four different sections. The first section consists of a listing of the setup program specifying the parameters to be used in the experimental and analytical part of the test. This is followed by the plot of enthalpy versus temperature recorded during the run The main exothermic reaction of interest occurs between 220'C and 440°C. The program sets these temperatures as limits of integration to quantify the scale of this exothermic peak. The DSC processor reproduces this portion of the curve as shown. Finally it prints the results of the integration and the kinetic ‘ parameters based on the regression analysis performed. The chief derived parameters are the following: order of reaction, n: Activation Energy, EA: frequency factor, f; heat of combustion, Q. -102- 54542_D§Q_29212!§!I The three main components of the equipment are : 1)The DSC measuring Cell: Mettler DSC20 2)DSC TA 3000 Processor 3)Swiss Matrix Printer £4§42;1_D§£_2211 The DSC Cell is outlined in Figure 4.7 showing the the two cells, probes and the relevant temperatures and heat flows. A111112.11HRIBLIEB£_QQEIBQL In heat flow calorimetry it is necessary that the reference temperature Tr follows a predetermined linear temperature program with heating rate T'. In reality this is possible only with an empty reference pan since a sample will likely show first-order transitions . In the case of a dynamic temperature program, the reference temperature necessarily lags behind the furnace temperature Tc as no heat flow, Q'r, to the reference side can take place without a temperature difference Tc-Tr. The lag is compensated in the TA3000 by means of a temperature advance as is shown in Figure 4.8. The measurement is terminated when the furnace temperature, Tc, reaches the preset final temperature. The DSC curve is only recorded -103- K I N C I I c A N A L V ‘ I 3 o-VIO-O. 11a.) ICAN Phaanttlnt 87.“! Y‘HP. .C 2.. ant! l/Hlu 20 (nu "fl'. ‘C 53. '2"! ISO. BIN. O PLOV CH 1. "Audi F. a“ 33 crest! “/0 3. Pan tvrt 118 1 LINK. -H . KIN('IC ANOLVfll' OVN/ISO 1i. 1 37A". 33. (no 4.. .fltILIH‘ 71" . ALPHA .700! 0.! ALPHA ‘HO 0.. PLO! en 1. rLot "001 I0! $INOL¢ “LPHR . GPPLXKD IINCVIC. 1.0/00! 1’: O IO‘HI. N0. llnvtnntunt 0c ”(It FLuu lxovutnnnL--) ’n. fl L; no». a I 1 DOI-.‘ 1 \ ”..- \R— fix t“. A J \ Isnrtunluut '6 "(0| ftOu Ix°nu£nnng~-) .. ... O. J .u (no -4 33:3,, .n 410 «01.n- Ptnu IIHP. 0c 338.) "Chet. can!“ 0.8. cont. glut? 0.0: t n IJIHOL 93.33 cant. slant a... L" K. $8.43 CONr. Llnlt 0.7. -------- "(Ytgtn In)... svsntn ~-------- -104- 00 Reference Pm O O DSC Samar Probe . ° 0 Hear C011 0 O O O o 0 Body of Furnace o O O O O O O O O O O O O O 0 -0 Tc Signal AT Skull 4—sze Gas Inlet Fig. 4.7 05: Cell. Tempcnfun T /T, - Db'phyed Tempennne Final Tempcnme ------------------- r. . "’ T, - Seraph Temperature , ‘ : Melting Pom! 71 lung]. met: I Tenpeniue unetpmn- ‘h‘ -T 47‘ T". E .' I ligaments-1 Om: Isomem um: r p e new 9590 phase Fig. 4.8 Temperature in the DSC measuring Cell as a Function of Time. -105- with respect to a reference temperature, Tr, which is lower by the amount of the advance. W The thermo-elements attached to the DSC sensor determine the temperature difference AT = T5 - Tr and produce a graph (see Figure 4.9). 11§1111_IEI_D§Q_EIQEBL The diagram of the test cell in Figure 4.7 shows that the heat flow to the sample , H', is equal to the difference between the two heat flows Q's and Q'r , where: Q's = Heat flow to sample pan, and Q'r - Heat flow to the reference. Thus, heat flow to the sample is equal to : H! g 9'5 - Q'r [4.2] According to the thermal analogue of Ohm's Law Q'= (T2 - T1) / Rth (i.e. the heat flow is proportional to the driving’ force, AT, and inversely 'proportional to the thermal resistance Rth)- Applying this concept to the DSC cell gives: -106- Dynamic Isothermal Isomennal Measun‘nz Measuring Phase Phase Phase I I 0 ' : > t, '1'; : I AT - T, - “r, I I AT I I I I w | l Fig. 4.9 AT Signal (Primary Signal) as a Function of Time, t, and the Reference Temperature, Tr. -107- H' = Q's ’ Q'r = Tc'Ts ‘ Tc'Tr [4-3] Rth Rth For reasons of symmetry, Tc and Rth have the same value for the sample and the reference. Thus: at = TS - Tr [4.4] The temperature difference , AT = Ts - Tr, is measured by the sensor thermocouples. From the thermocouple equation, AU = T.S, it follows that: H' = AU / Rth . s [4.5] The two terms in the denominator of Equation [4.5] are functions of the actual temperature and can be combined to define the calorimetric sensitivity : E = nth . s [4.6] B may be divided into a temperature dependent (relative) term Erel and a temperature independent term EIn. Specific to the measuring cell : -108- E = Brel - EIn [4.7] Thus the heat flow to the sample is equal to: H' = AU / EIn - Erel [4-8] The temperature dependence of Erel is contained in the TA processor as a polynomial : Erel = A + ST + CT2 [4.9] The specific sensor parameters A, B and C are fixed with the coding plug in the measuring cell. EIn is determined by calibration using the known heat of fusion of Indium. EIn corresponds to the coding plug setting "Medium Sensitivity" with the standard sensor (c. 11 uV/mW). The ‘ primary signal is converted once per second using Equation [4.46] for the on-line plot of the printer / plotter. No account is taken in Equation [4.46] of the fact that the heat capacities as well as the resistances are dependent on the path of the heat flow. They cause damping with a specific time constant similar to an electrical RC term. The time constant, tsignal (c.7.5s), leads to a broadening of the DSC signal. The original heat flow to the sample may be reconstructed with the help of the following equation ( deconvolution ): -109- H. = At] + tsignal 0 U. [4.10] EIn - Erel Theoriginal heat flow is calculated in all off-line evaluations, printed on the printer/plotter and is used for partial and total integrals. 1L§&1&§_HEIILEB_£BQQE§§QB The function of the Mettler processor is: 1) To enter information necessary for the operation of the DSC measuring cell: 2) To control the furnace in the measuring cell; 3) To acquire and store the curve data for various evaluations; 4) To analyse the measured curve using various evaluation methods and calculate the final numerical results; and 5) To provide an interface for the printer to print out the experimental parameters, measured curves, calculated curves and results. -110- 1;§;11£_£BIEI§B_L_£LQII§B The printer/plotter is interfaced with the TA processor, allowing for the retrieval of information (calculations, graphs etc.) from the processor memory. Aii11_SLLI§BAIIQE_QI_IZ§_2§£ A series of parameters are specified in the list of configuration data. They are needed for the measurement and control of the furnace temperature, and for converting the T signal to the heat flow H'. The relevant values for the measuring cell and the sensor are automatically read when the measuring cell is connected. The measured data is divided into accessible and protected data. The accessible data are displayed with the CONFIG function if they are to be altered during calibration. The calibration of the DSC is divided into : 1) Heat flow calibration 2) Temperature calibration. The calibrations are carried out using two different procedures. -111- 1) Heat flow calibration The heat flow is calibrated using the heat of fusion of a known quantity of Indium. The calorimetric sensitivity, EIndiumr is subsequently entered in the list of configuration data. The BIndium value is determined as the average of a number of calibration runs. 2) Temperature calibration A standard pan is supplied with the DSC for calibration. The pan consists of’ a known quantity of indium, lead and zinc in separate compartments. The coefficients of the sensor temperature-dependence equation (A, B and C) are obtained from the fusion of the three metals, and are entered automatically. Calibration data are entered into the configuration data both manually and automatically. Wm The time constant for the temperature equilibrium between the furnace and. the DSC‘ sensor, tlagr are determined experimentally. For this purpose, the melting point of indium is determined at different heating rates using the -112- purity method. Thus, for A > B, where TA and T3 are the melting points at heating rates A and B K/min, respectively. Table 4.3 shows the results of the calibration exercise: 1§§13_1;; DSC Calibration data. EIAIIH§_BAIB 1_szin M.P. Indium (1) 174.2 175.3 M.P. Indium (2) 174.6 175.3 MEAN 174.4 175.3 Thus tlag = 22 + 60(175.3 - 174.4)/10 = 27.4 secs This value is entered into the configuration data. Aiiii_flAHILE_2BEELBAIIQE_LED_IEEEBIIQE The standard operational procedure outlined in the Mettler handbook is used to prepare and insert the sample. The thermal data for the run are entered via the keypad: the -113- instructions for the curve or output data analysis are also be entered. W The DSC may be used to analyse substances under the following headings: - Screen - Kinetic - Purity - Integ - Cp, Specific Heat Different programming sequences are required to implement the respective temperature programs and perform the attendant analyses. Details of these programming procedures are listed in Appendix II. -114- OR ' LO KENT O - £11.131392991193 As outlined in Chapter 3, much of the research conducted to date on the question of self-ignition of powders is empirical. While some attempts have been made to propose theoretical solutions (cf. the Frank-Kamenetskii, Semenov models ), the basic heat transfer assumptions used tend to oversimplify the problem. Such simplification was needed, however, due to the difficulty in handling the mathematical solution process. As Figure 3.2 shows, the real situation entails finite surface and internal resistance to heat transfer, as well as internal heat generation. By using a finite element technique to solve this problem, an elemental approach is possible to model the heat transfer conditions. A numerical analysis can then be used to simulate the heat transfer over time. Finite differences or finite elements (F.E.) may be used to model the heat transfer mechanism of interest here. F.E. was the preferred method as it facilitates the extension of the model to a range of geometries including amorphous shapes. The F.E. model is first used to simulate the surface and internal heat transfer phenomena. After '115- verifying this model, the heat generation term is added. Further verification requires experimental results. The equation governing the heat transfer / temperature profile in a material which is undergoing a self-heating reaction can be stated as ( Carslaw and Jaeger, 1959 ) : pc aT/at = k 0217:”:2 + G [5.1] where p - material density ( kg / m3 ) O I material specific heat capacity ( J / kg K ) - thermal conductivity ( W / m K ) temperature ( K ) - time ( s ) Q r? *3 K I - heat generation term ( W / m3 ) In the case of a degradation / oxidation reaction the heat generation term, G (W/m3), may be expressed in terms of the following general equation (Drysdale, 1985): G = Q cin r e('EA / R T ) [5.2] where Q - heat of combustion ( J / mole ) -116- n = order of reaction Ci = concentration ( mole / m3 ) f = a pre-exponential factor, units depending on order of reaction, n EA - Activation Energy ( J / kg mole ) R a Gas constant, 8.3143 J / K mole T = temperature ( K ) The initial temperature is assumed uniform in the sample : T = Tinitv t = 0 [5.3] To allow for the widest range of environmental conditions, a convective boundary condition is specified : k aT/an = h ( T - Tamb ) [5.4] where n is the vector normal to the surface of the reacting mass, Figure 5.1. We}: The first step in developing the model is to decide upon a solution technique for Equation [5.1] which allows variation of as many experimental conditions as possible. The Finite ZElement approach. allows ‘variation, both in terms of system geometry and composition. The technique '117- Fig. 5.1 n is the Vector Normal to the Sphere Surface. -118- also permits a fine grid to be constructed, in space dimensions and time steps. Thus, many nodes / elements can be monitored near a sample surface where the critical heat loss to the environment takes place. Alternatively, attention may be focussed on the system hot-spot i.e. sample center or near-center. As the solution technique is iterative, the time increment can be altered to allow detection of the 'take-off point' at which the sample temperature increases very rapidly as the runaway reaction commences . Anderson et al.( 1974 ) proposed a solution technique for a simplified version of equation [5.1] using the Finite Element approach. The authors sought to predict values of the Frank-Kamenetskii parameter, 5cr: using theoretical temperature and property data. Section 3.3.2 contains the details of this study. R DU Equation [5.1] without the heat generation term, i.e. G = 0. This allows verification of the heat transfer part of the model, involving the heat conduction, convection and capacitance terms. The simpler model has application in terms of the tests used to establish the environmental conditions in the oven. -119- Spheres of sample material are typically used in a laboratory test oven to establish the commonly used quantities of MIT & 5cr- The symmetry allows the mathematics to be reduced to a uni-dimensional problem, and thus the study can concentrate on the key factors which influence the start of ignition viz. product composition, ambient temperature and critical sample dimension. Given the nature of the Finite Element technique, subsequent alteration of the shape functions allows considerable generalisation of the model and its applications. In contrast to the Semenov or Frank-Kamenetskii models, the present study includes no simplifying assumptions in terms of either surface or internal resistance to heat transfer. Both thermal resistances are assumed to have a ' finite value i.e. Figure 3.28. MW With the problem reduced to the axi-symmetric, uni- dimensional, no-heat-generation problem, Equation [5.1] reduces , in spherical co-ordinates, to the following (Carslaw et al., 1959) : -120- kr a2 T + 2 kr a T = pc 0 T [5.5] or2 r a r a t The boundary equation is equal to Equation [5.4] with n = r, the radial length co-ordinate; and kr is the radial thermal conductivity. The initial condition is given by Equation [5.3]. Thus, the problem is to calculate the temperature profile over time of a sphere of radius r, initially at a temperature Tinitr exposed to a temperature of Tamb- The solution of the problem, using the Finite Element technique is outlined below. Thus, the time / temperature profile and the combustibility parameters can be accurately estimated. iizi1iz_zinife_zlenent_§:i§ The sphere is first divided into a number of connective elements along a typical radius as shown in Figure 5.2. 1 2 3 i j E-l E E+1 119321.512 Finite Element Grid. Additional elements are included near the sphere boundary as the surface heat transfer dictates the rate of heat loss from the sample. It is the heat loss which dissipates -121- the heat generated by the oxidation reaction. Hence, particular attention needs to be paid to the temperature profile near the boundary. The solution is obtained by finding a function T (r,t) , satisfying the boundary conditions , which minimises an integral quantity called a functional. The functional incorporates both the field equation [5.5] and the boundary condition [5.4] , is minimised with respect to time to give the approximate temperature profile of the sphere. The integral term is I = 1/2 [ kr(T')2 + 2ch aT/ot]dV V X Y + 1/2 h(T-Tamb)2ds [5.5] S 2 . where T' = aT/ar and the integrals are the volume and surface integrals, respectively [ Segerlind, 1976 ]. As the sphere is divided into E elements, the integral can be evaluated separately for each element 'e' and the results summated [ Myers, 1971 ]. Thus : I a 1(1)+I(2)+I(3)+...I(e)...+I(E-1)+I(e) [5.7] Taking into account the three modes of heat transfer / -122- loss in a typical element, the integral 1(9) may itself be subdivided into conductive, capacitive and convective sub- integrals as follows: I(e) = Ik(e) + Ic(e) + Ih(e) [5,3] The subintegrals correspond, respectively, with the expressions labelled as X, Y and z in Equation [5.6] above. When each of these integrations is performed and the integral is minimised with respect to the element temperature T(e), the 'element equations' describing the heat transfer result. The details of the calculus operations are included in Appendix I; only the results are included below. WW [K(e)] = 4 t k(e)(rj3-rj3) [ 1 -l] 3 (rj-ri)2 [5.9] iiz11iA_I1enent_saeesiienee_eetrix [013(3)] ”Lt—M):- C11 c12] 60 (ri-rj) c21 c225 [5.10] where 611 = 2rj5-20rj2ri3+30rjri4-12r15 -123- 612 = 3rj5-5rj4ri+5rjri4-ris C21 = C12 c2; erj5-30rj4ri+20rj3r12-2r15 WW3: This matrix is null for all elements except for the Eth or last element where it has the form [HE] = 4 x h o o o R2 [5.12] Summarizing the element matrices for each heat transfer component gives the global conduction, capacitance and convection matrices. Finally, the system force ‘vector is computed. In the absence of a heat generation term, this vector consists solely of the convective force at the surface and may be written as: (F) = 4 t h Tamb [5.13] N0'°OOO -124- Thus, (F) is a column matrix with one non-zero term. WW When the element matrices are assembled for all elements, the equation may be written in global matrix form: {K} {T} + {C} {T}' = {F} [5-14] where (K) incorporates the conduction. as well as the convection element. matrices, and (C) includes all the capacitance terms. (T)~ is the column matrix of temperature/time derivatives ( i.e. a typical row being oT(e)/at). Solution of equation [5.14] gives the time / temperature profile for the conditions outlined. W The Crank-Nicolson technique (Crank et al., 1947) , was used to solve Equation [5.13], approximating’ the time differential over a time interval At. This gives rise to the following system of matrices : [{K}+(2/At){C)]{Tnew} = [(Z/At){C)'{K}]{Told} + 2{F} [5.15] -125- As the derivatives are evaluated at the midpoints of the time interval, the nodal temperatures can be similarly determined (Segerlind, 1976). A similar approach was used by zienkiewicz et al. (1970) in solving the transient field problem using finite elements. Taking this into account, Equation [5.9] may be written in a compact form suitable for an iteration process: {A} {T)new = {B} {T)old + {F} [5-15] where (A) - (K) + {B} and (B) - (2/At) {C} The resultant {T}new is the temperature profile after a time of At/Z. iiZ12_IIBIIIQAIIQE_QI_IEI_E£LI_IBLEEIIB_HQD£L i121211_QQflBLBLEL£_BELLIII£AL_HQDBL In order' to 'verify‘ the accurate. working’ of the heat transfer model developed above, an analytical solution of a test problem was used as a base line. A program was written to evaluate the center-point temperature profile for a sphere with finite surface and internal resistances to heat transfer. This was based on the following equation, Wong (1977): -126- T - Tamb - 4 sin Bn-Bncos an exp(-[Bn/R]2) a t Tinit'Tamb ZBn' SinZBn [5.17] where 8n are roots of En cot Bn = (1-hR/k) R - sphere radius a a thermal diffusivity £121212_§LHELB_£EQ§L£H To compare the two solution methods , a set of data was assembled for an aluminium sphere. The data were for an aluminium sphere used to establish test parameters in one of the experimental ovens. ALHMIHIQM_§£EEBE Radius 0.0508 m Specific heat capacity 900 J/kg K Thermal conductivity 220 W/m K Density 2730 kg/m3 IE§I_§QHDIIIQH§ Surface convection coefficient 13.75 W/m2 Tinit 298 K Tamb 423 K ii2i211_IEHIEBLI!BI_RBQIILE§ The results of the analytical and numerical solutions of the temperature profile of the sphere are presented in Table 5.1 and plotted in Fig. 5.3. Fig. 5.4 shows an expanded view of part of the plot. -127- 13§L£_§;1 Analytical and numerical solutions for the aluminium sphere. -128- iii—21m M231 211mm (mins) Temp.(’C) Temp.(°C) 0 25.00 25.00 1 27.31 27.35 2 29.69 29.76 3 32.03 32.12 4 34.32 34.43 5 36.57 36.70 6 38.77 38.92 7 40.94 41.10 8 43.05 43.23 9 45.13 45.33 10 47.17 47.38 11 49.17 49.39 12 51.13 51.37 13 53.05 53.30 14 54.93 55.20 15 56.78 57.00 16 58.59 58.88 17 60.37 60.67 18 62.11 62.42 19 63.82 64.14 20 65.49 65.82 21 67.13 67.47 22 68.74 69.09 23 70.32 70.68 24 71.87 72.24 25 73.39 73.76 26 74.88 75.26 27 76.34 76.72 28 77.77 78.16 29 79.17 79.57 30 80.55 80.95 31 81.90 82.31 32 83.22 83.63 33 84.52 84.93 34 85.79 86.21 35 87.04 87.46 36 88.26 88.69 37 89.46 89.89 38 90.64 91.07 39 91.79 92.23 40 92.92 93.36 41 94.03 94.47 42 95.12 95.56 43 96.19 96.63 Tine Zinite_lleeent (mins) Temp.('C) Temp.('C) 44 97.23 97.68 45 98.26 98.70 46 99.26 99.71 47 100.25 100.72 48 101.22 101.70 49 102.16 102.60 50 103.09 103.50 51 104.00 104.50 52 104.90 105.30 53 105.78 106.20 54 106.63 107.10 55 107.48 107.90 56 108.30 108.70 57 109.11 109.60 58 109.91 110.30 59 110.69 111.10 60 111.45 111.90 61 112.20 112.60 62 112.93 113.4 63 113.65 114.10 64 114.36 114.80 65 115.05 115.50 66 115.73 116.20 67 116.40 116.80 68 117.05 117.50 69 117.69 118.10 70 118.32 118.70 71 118.94 119.40 72 119.54 120.00 73 120.13 120.50 74 120.71 121.10 75 121.28 121.70 76 121.84 122.20 77 122.39 122.80 78 122.92 123.30 79 123.45 123.80 80 123.96 124.40 81 124.47 124.90 82 124.97 125.40 83 125.45 125.80 84 125.93 126.30 85 126.40 126.80 86 126.86 127.20 87 127.31 127.70 -129- ('C) Temp. Fig. 140 120 100 l A 1 L l 50 100 150 Time (mins) 5.3 Analytic vs Finite Element Solution for Aluminium Sphere Heat-up. -130- 140 " Finite Element Solution ('0) Analytical Solution 130 " Temp. 125' J A J 10 m 90 100 110 120 Time (mine) Fig. 5.4 Analytic vs Finite Element Solution for Aluminium Sphere Heat-up, (Expanded View). 120 < . -131- For the Finite Element program a grid of fourteen elements is used, with more and smaller elements near the sphere surface. A time step of 120 sec proved to be satisfactory. No stability problems are encountered using the Gauss- Seidel solution routine with Equation [5.14] when calculating the nodal temperatures at each time step. Using the analytical technique to calculate the center- point temperature as a function of time and conditions, the series solution, Equation [5.17], is limited to five terms. The relevant 8n roots are : 0.09755, 4.4941, 7.7257, 10.9044, 14.064, (Abramowitz et al., 1972). The fact that the graphs for F.E. and the analytic solution practically coincide , shows that the level of accuracy attainable with the Finite Element program is satisfactory._ Table 5.2 and Fig. 5.5 show an analytic and numerical comparison for a different problem. The sample is an apple (Data: k - 0.418 W/m K, c - 3.766 kJ/kg K, p - 787 kg/m3, R - 0.04 m, h - 60 W/ m2 K) being cooled from 30'C by a O'C airstream. Again the Finite Element approach yields very good agreement with the analytical solution. Varying the surface heat transfer conditions can be modelled at will with the F.E. program as can be seen in Fig. 5.6 .and '132- Ta21g_§;z Analytical and numerical solution for a cooling problem. BELLXIIELL ZIEIIE_£L§§EEI 1133 IEfiEEBAIEBI 1133. IEEEEBAIEBB (mins) (’C) (mins) ('C) O 30.00 0.5 30.00 1 30.19 1.0 30.00 2 30.00 1.5 30.00 3 30.00 2.0 30.00 4 30.00 2.5 30.00 5 29.99 3.0 30.00 6 29.96 3.5 30.00 7 29.88 4.0 30.00 8 29.72 4.5 30.00 9 29.47 5.0 30.00 10 29.11 5.5 29.99 11 28.66 6.0 29.97 12 28.11 6.5 29.94 13 27.47 7.0 29.90 14 26.78 7.5 29.85 15 26.02 8.0 29.77 16 25.23 8.5 29.67 17 24.41 9.0 29.55 18 23.58 9.5 29.41 19 22.73 10.0 29.25 20 21.89 10.5 29.06 21 21.05 11.0 28.85 22 20.22 11.5 28.61 23 19.41 12.0 28.36 24 18.62 12.5 28.08 25 17.84 13.0 27.79 26 17.09 13.5 27.47 27 16.35 14.0 27.14 28 15.66 14.5 26.80 29 14.98 15.0 26.44 30 14.32 15.5 26.08 31 13.69 16.0 25.07 32 13.08 16.5 25.31 33 12.50 17.0 24.92 34 11.92 17.5 24.52 35 11.41 18.0 24.12 36 10.90 18.5 23.71 37 10.41 19.0 23.30 38 9.94 19.5 22.89 39 9.48 20.0 22.48 40 9.06 20.5 22.06 41 8.64 21.0 21.65 42 8.25 21.5 21.24 -133- W W 11!! 15525381233 1153. 15525351233 (mins) (‘C) (mins) ('C) 43 7.88 22.0 20.84 44 7.52 22.5 20.43 45 7.18 23.0 23.03 46 6.85 23.5 19.63 47 6.54 24.0 19.24 48 6.24 24.5 18.85 49 5.95 25.0 18.46 50 5.68 25.5 18.08 51 5.42 26.0 17.70 52 5.17 26.5 17.33 53 4.94 27.0 16.97 54 4.71 27.5 16.61 55 4.50 28.0 16.26 56 4.29 28.5 15.91 57 4.09 29.0 15.57 58 3.91 29.5 15.23 59 3.73 30.0 14.90 60 3.56 30.5 14.58 61 3.39 31.0 14.26 62 3.34 31.5 13.94 63 3.09 32.0 13.64 64 2.95 32.5 13.34 65 2.81 33.0 13.04 66 2.69 33.5 12.75 67 2.56 34.0 12.47 68 2.45 34.5 12.19 69 2.33 35.0 11.92 70 2.23 35.5 11.65 71 2.12 36.0 11.39 72 2.03 36.5 11.14 73 1.93 37.0 10.89 74 1.85 37.5 10.64 75 1.76 38.0 10.40 76 1.68 38.5 10.17 77 1.60 39.0 9.94 78 1.53 39.5 9.71 79 1.46 40.0 9.49 80 1.39 40.5 9.28 81 1.33 41.0 9.07 82 1.27 41.5 8.86 83 1.21 42.0 8.66 84 1.16 42.5 8.47 85 1.10 43.0 8.27 86 1.05 43.5 8.09 '134- “F .. 30 "\ .0 a. "3 E 20 - h a 1‘3. ‘4 10 - o l o 50 100 150 Time (mins) Fig. 5.5 Analytic vs Finite Element Solution for Apple Cooling. '135- 13p1g_§;1 Time/temperature data for varying convection conditions. 1189 (mins) omflmmbUNi-J h=100 Slim. (’0 30.00 30.00 30.00 30.01 30.02 30.04 30.05 30.05 30.04 30.00 29.94 29.84 29.72 29.57 29.38 29.17 28.94 28.68 28.40 28.10 27.79 27.46 27.12 26.77 26.41 26.05 25.68 25.31 24.93 24.56 24.18 23.80 23.43 23.05 22.68 22.31 21.94 21.58 21.22 20.86 20.51 20.16 19.82 h=200 1289... ('C) 30.00 30.00 30.01 30.02 30.04 30.07 30.09 30.09 30.06 29.99 29.88 29.72 29.51 29.25 28.94 28.60 28.22 27.80 27.36 26.90 26.41 25.91 25.40 24.88 24.35 23.82 23.29 22.76 22.23 21.70 21.18 20.67 20.16 19.65 19.16 18.68 18.20 17.73 17.27 16.82 16.38 15.95 15.54 -136- h=600 me... ('C) 30.00 30.00 30.01 30.04 30.08 30.13 30.16 30.15 30.08 29.94 29.71 29.39 29.00 28.53 28.00 27.42 26.78 26.11 25.42 24.70 23.97 23.23 22.49 21.75 21.01 20.29 19.58 18.88 18.20 17.53 16.88 16.25 15.64 15.04 14.47 13.91 13.37 12.85 12.35 11.87 11.40 10.96 10.52 Time (mins) 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 h=100 Temp, (° C) 19.48 19.15 18.82 18.49 18.17 17.86 17.54 17.24 16.94 16.64 16.35 16.06 15.78 15.50 15.23 h=200 Em; ('C) 15.12 14.72 14.33 13.95 13.58 13.22 12.86 12.52 12.18 11.85 11.53 11.22 10.92 10.62 10.34 '137- h=600 ('C) 10.11 9.71 9.32 8.95 8.60 8.26 7.93 7.61 7.31 7.02 6.74 6.47 6.21 5.96 5.72 40 " A 30 “"24. p ._‘..'.. 3' e Z)” l- 10 ’ 0 - l 4 l + l 4 l A 0 Z) 40 El 80 1(1) 120 140 Time (mins) Fig. 5.6 F.E. Technique used to Model Diiierent Values of Heat Transfer Coefficient. '138- Table 5.3. As the data shows, the model produces satisfactory results under varying conditions. The next step in the solution procedure is to include the Arrhenius heat generation term in the global statement of the problem and to solve the resulting equation. U 0 T The most general expression of the heat generation term is given in Equation [5.2]. If the Frank-Kamenetskii reactant assumption is accepted (i.e. the rate of this exothermic reaction is independent of concentration of reactant) the generation term may be simplified to G = Q f e( '31 / R T ) [5.18] where Q 8 heat of combustion ( J/m3 ) f - frequency factor Equation [5.18] differs from Equation [5.2] only by the fact that the term Ci” has an index of zero i.e. -139- The model is developed including a variable value for the index to allow the reaction kinetics to be properly accounted for in the simulation exercise. The case of a zero concentration index, or zero-order reaction, will be studied as a special case to permit a comparison with the more traditional models used in this area. Anderson et al. (1974) proposed a method of solving the overall Equation [5.1] for a zero order reaction. It has been found experimentally that during the course of an oven test, significant consumption of reactant occurs. The rate of the consumption influences Equation [5.2], the generation term, and hence influences the transient solution of the overall heat transfer Equation [5.1], and its ability to properly identify parameters such as the induction time and the MIT. The model developed here should account for the non-zero reactant condition. The temperature profile which results when the heat generation term is included, may cause instability problems when the conditions cause a runaway reaction. With the onset of thermal runaway, the temperature increases very rapidly approaching the point of ignition. Thus, the term (T)', the time derivative of the temperature, will increase very rapidly, causing the numerical computation to become unstable. Care must be taken with the size of the time interval of the -140- differentiation when a significant temperature rise is detected, with smaller and smaller time steps necessary as the MIT point approaches. i12121Z_IEI_ZQB92.!IQIQB_IL§EEEI_HBIBIX In the finite element formulation, the integral term as given by Equation [5.6] is rewritten to include the heat generation term, giving : I = %[kr(T')2 + 2ch aT/at - ch ] dV v + Isl/2 h [T-Tamb]2 as [5.19] Extracting the heat generation term from Equation [5.19], G = Q cin f e'EA/RT [5.2] gives for element (e) : 13(8) - Ive Tie) dV = 4 T rj G T(e) r2 dr ri ( Note : G is evaluated at the end of each time interval, as a function of r and the nodal temperatures Ti and Tj.) '141- The integral then may be written as : 16(9) = 4 1 rj r2 G [NiTi + NjTj]dr ri 15(9) is differentiated with respect to T(e) to minimise the component of the 'functional'. The result yields the two components of the generation term : 91(9) = 4 x rj r2 G Ni dr ri 4 1 U rj e‘EA/RT(r2rj-r3) dr rj-ri ri ( 91(9) is the contribution of element e to the temperature at node i) and gj(e) = 4 t J rj r2 G Nj dr ri = 4 x U I rj e’EA/RT(r3-r2ri) dr rj-ri ri [5-201 ( 93(9) is the contribution of element e to the temperature at node j ) where U = Q f, with Q and f defined in Equation [5.18]. '142- As no analytical solution is possible for the integral term due to the presence of the Arrhenius exponential term, e’EA/RT, it is evaluated after each step of the iteration process. In the program listed in Appendix III, this operation is performed using a 32-point Gaussian quadrature procedure. The elemental terms 91(9) and gj(°) combine to give the heat generation component contribution to the Force Vector, { G ), in the following manner: 1 91 921+922 ( G l = . [5.21] 98E’1+QEE L 9E+1E After each time increment, the ( G } vector is calculated and included in the global Equation [5.15] . As ( G } must be calculated anew after each time iteration, it is convenient to include the heat generation term as part of the force vector. This represents the volume integral of the gT term in Equation [5.15]. The global equation is solved in the manner previously outlined. -143- As outlined in Section 4.2.1, the surface heat transfer coefficient is calculated from the time/temperature profile of an aluminium sphere for each of the sample radii used in the experimental phase. The standard equation is used: T - Tamb - e-(3h/Rpclt Tinit - Tamb [6.1] The sphere data used in the equation for a typical 10.16 cm (4") diameter sphere are: Ambient temperature : 143.2 'C Initial temperature : 25 ‘C Radius (R) : 0.0508 m Density (p) : 2730 kg/m3 Specific Heat (C) : 900 J/kg 'C Table 6.1 shows the data for the time/temperature profile which is shown plotted in Figure 6.1. Equation [6.1] is used to evaluate the effective h value at each of the recorded data points. The results of these calculations are also shown in Table 6.1. The outlying points are eliminated from this data set. These are chiefly associated with the start-up period. The average h-value and standard deviation are then evaluated. -144- for 4” 11.111193221331112! (mins) 1.20 1.76 2.26 3.31 4.89 6.21 7.13 8.04 9.11 10.13 11.51 12.47 13.24 14.25 15.58 16.72 17.82 18.99 20.26 21.92 23.06 24.49 25.91 27.41 29.07 30.90 32.63 34.67 36.83 38.89 41.08 43.07 44.65 46.89 49.23 50.98 ('C) 29.41 31.25 33.11 36.60 41.26 45.03 47.72 50.42 53.01 55.69 58.48 61.16 63.83 65.91 68.11 70.67 72.76 75.23 77.11 80.06 82.42 84.96 87.13 89.19 91.82 93.89 96.41 98.94 101.37 103.98 106.04 108.10 109.58 111.25 113.12 114.04 -145- 133L§_§;1 Time/temperature and 'h'-value data sphere. Obi-zllng (W/MK) 21.93 21.36 21.77 21.60 20.94 20.69 20.72 20.84 20.55 20.54 20.02 20.27 20.82 20.63 20.15 20.21 20.10 20.16 19.86 19.80 19.96 20.00 19.92 19.77 19.83 19.58 19.65 19.61 19.52 19.63 19.49 19.51 19.49 19.31 19.24 19.00 ('0) Temp. 12°F b . . I. I 100 '- . ' .I h .I I .I 80 " .e 'I h .- I I 60 - _' b ... 4o .- 20 l A j l A l l A I O 10 20 30 4O 50 60 Time (mine) Fig. 6.1 4" Aluminium Sphere Heat-up Curve. -146- Thus, for the nominal 4" sphere, a value of 20.03 W/mZK was determined, with a standard deviation of 0.48 over 31 points on the curve shown in Figure 6.1. Values of the surface heat transfer coefficient were similarly obtained for 3", 2" and 1.5" spheres and their respective values are recorded in Table 6.2. 15212_§&z Experimentally determined values of the surface heat transfer coefficient W _'hi (inches) (W/mZK) 4 20.03 3 20.12 2 21.26 1.5 21.78 As the values in Table 6.2 indicate, the fixed fan speed, corresponding to a motor voltage of 210V, gives, on the whole, a consistent value of heat transfer at the sphere surface. A slightly increased value for 'h' was found with decreasing sphere size. These figures represent the effective surface heat transfer coefficients and are used in the simulation exercise to plot the time/temperature profiles for the respective experimental sphere sizes. The -147- good correlation between the experimental and simulated profiles, as may be seen, for example, in Figures 6.5 - 6.8 below, proves the quality of these values. i11I1_HII_DAIA_IQBAAEQEHQBB_§EIH_HILK_RQ!DEB The relevant experimental centerpoint time/temperature plots for the 4" sphere are shown in Figure 4.5. The lowest ignition and highest non-ignition temperatures give rise to the critical centerpoint temperature profiles shown as Figures 6.2, 6.3 and 6.4 for the 3", 2" and 1.5" samples of skim powder, respectively. Table 6.3 shows the ambient temperatures and the resultant 'averaged' Minimum Ignition Temperature for each sphere size for the Avonmore skim milk used in the tests. In the table, the Ignition column is the lowest ambient temperature at which the sample ignited. The figure for No Ignition is the highest ambient temperature at which the sample failed to ignite. The difference between these two figures is listed in the 'Accurecy' column. Thus in the case of the 4" sphere, there is a span of 1.63‘C between ignition and non-ignition. This rises to above 4'C for the 3" sphere. This is quite common with the oven test technique, where there may be a gap of some 10'C between -148- 8m) ('0) 4m) Temp. 2m) Fl I 146'C 's i I i i - i [i .; 141.8'C .II L. I I I I A l A l 1000 2000 Time (mine) Fig. 6.2 3” Sphere Milk Powder Ignition Curve. -149- Temp. ('C) 162 'C 160.4'0 I L L A l A J Z!) 400 6m 8m 1m Time (mins) Fig. 6.3 2" Sphere Milk Powder Ignition Curve. -150- 1000 ('0) Temp. 173.8'C 171.5'C 300 400 500 Time (mins) Fig. 6.4 1.5" Sphere Milk Powder Ignition Curve. -151- non-ignition and ignition in some instances (Synnott et al., 1986). 11p11_§L1 Radius/Experimental MIT data for Avonmore skim milk. 0 19311198 82.19811128 811 AQEHIQEI (Nominal) ("l ('C) ('C) ('C) ('C) 4 138.90 137.27 138.08 1.63 3 146.00 141.82 143.91 4.18 2 162.00 160.37 161.19 1.63 1.5 173.82 171.46 172.64 2.36 In using the oven technique to assess products for criticality, it is important to gauge the possible effect of experimental method on the result. temperature is the key parameter of interest. In the experimental set-up this is measured using thermocouples. Larkin (1984) proposed a formula to account for the heat lost through the thermocouple wire when measuring thermal diffusivity in a cylinder/can. of product. Larkin's equation effectively calculates the correction factor due to the thermocouple. Thus : T'(t) = T(t) (1 + 2e'mR(tanh(mR)-1)/(mR)2 +2(l/mR - tanh(mR))/mR)) [6.2] For ‘the [experimental conditions of interest. here, the above variables take on the following values: -152- T'(t) = adjusted temperature at center of sphere at time t T(t) = measured sphere center temperature at time t m - (hP/kA) 1/ 2 where h - 19.21 W/mzK (convection coefficient) P - 6.28 X 10’4 (perimeter of probe) k - 46.28 W/mK (probe thermal conductivity) A - 3.14 X 10’8 (cross-sectional area) R = 5.05 x 10'2 m (sphere radius) Calculating as per Equation [6.2] yields the result: T'(t) = 0.98 x T(t) Thus a maximum error of 2% is due to the thermocouple. As the aim of the measurement is to identify the ignition point, i.e. a. very large increase in temperature, this error is not significant. LAW To determine the kinetic parameters for the spontaneous ignition , a Bowes plot of 1n (6c, Tcrz/roz) vs. 1/Tcr was constructed based on the data listed in Table 6.3 -153- using a value of 3.32 for 5cr( Drysdale, 1985) , as outlined in Chapter 3, Figure 3.1 ( see Figure 6.5) . The coordinates for the plot are shown in Table 6.4. 13p11_g.1 Bowes plot data for Avonmore skim milk. :9 w; MgIgg-ZLIQZ-J. (cm) I/KTI 5.05 2.43x10'3 19.21 4.02 2.401(10'3 19.69 2.73 2.30x10'3 20.55 2.14 2.24x10“3 21.09 The data plotted in Figure 6.5 yield the following results: Slope :-9540 Intercept : 42.48 Correlation Coeff. : 0.995 Thus, EA/ R a 9540 , giving EA - 79.316 kJ/kg mole Similarly, ln(EAQf/kR) - 42.48 Thus, EA and R are known. The product thermal conductivity , k, for this sample can be found using the Maxwell-Euchen equation (MacCarthy, 1983), Equation [3.40]: Re 3 kair 1 ’ fV( 1 'biksol/kair} ) 1 + fv ( b - 1 ) [ 3.40 ] Firstly, product density must be calculated to derive a figure for fv, the solids fraction. -154- Ln CrTCr/r 01M 0.0023 0M4 1/Tcr ("0 Fig. 6.5 Bowes Plot. -155- £511211_DIE§IIX Milk powder density is determined from the compositional data by adding the component densities as outlined in Equation [3.52]. These components for the sample in question , together with their respective densities , are: gompgnent Content Density (It) (kg/m3) Moisture 4.25 1000 Fat 1.49 930 SNF 94.26 1600 Equation [3.52] gives a value for the powder density : Ppowder = 1544 kg / m3 The experimental bulk density value found for the powder was 600 kg/m3 (Section 4.2.2). This gives a porosity value of P = 1 - pbulk/psolid = 1- 600/1544 = 0.6114 [ 6.3 ] The 600 kg/m3 figure for bulk density is the figure used for' both the heat capacity term, in conjunction with Equation [5.10], and in estimating the heat generated per unit volume in Equation [5.2]. It is in line with the 610 kg/m3 figure quoted by Carr (1976) for milk powders. '156- To ‘use the ZMaxwell-Euchen expression for thermal conduction as outlined in Equation [3.40], the solids fraction must first be calculated. Given a porosity of 0.6114, fv is : The parameters 1J1 Equation [3.40] are thus assigned the following values : kair : 0.030 W/mK ksol : 0.419 W/mK fv : 0.3886 (Heldman et al.(1981)) From Equation [3.40], b 3 3 kair / ( 2 kair + ksol ) Thus: b = 0.188 This gives an effective thermal conductivity of: ke = 0.0716 W/mK '157- Returning to the calculations from Figure 6.5 yields a value for Qf, 'heat of combustion' = 2.11 X 1013 W/m3 Table 6.5 summarises the known and derived data for the Avonmore skim milk studied in the experimental phase as a reference to test the accuracy of the finite element model. Iablg_§&§ Experimental and derived milk-powder sample data. 231389.139: 5.1831921 231118 Specific heat Cp 1547 kJ/kg K Thermal Conductivity k 0.0716 W/mK Density p 706 kg/m3 Activation Energy EA 79.316 kJ/kg mole 'Heat of Combustion' Qf 1.1572 X 1013 W/m3 £12_D§£_DAZL Samples of skim milk powder were run in a Mettler Differential Scanning’ Calorimeter; By ‘using fixed. rate heating, the following parameters were determined for the powder : Heat of combustion, Q, in kJ/kg, the frequency factor, f, in s‘1 and the order of the reaction, n. These parameters are as previously defined in the basic heat generation equation, Equation [3.2]. They characterise .- the combustibility and nature of the -158- exothermic reaction under investigation. The data is used as the starting point for the simulation of the spontaneous ignition phenomenon thereby circumventing the time consuming empirical oven/hot-plate approaches to categorising powder combustibility. The procedure followed is: 1) Determine Q and f for a powder. 2) Simulate the heat transfer/generation process via the Finite Element model and predict the powder MIT and Time to Ignition. 3) Compare the MIT values from 2) with experimental oven results and with published data. é1Z11_QaLQBIHIIBIQ_DLIL_ZQB_HILE_RQ!DBB§ The results obtained using the kinetic analysis program on the DSC for skim milk are printed out per Figure 4.6. For the Avonmore skim milk the relevant DSC data are: AH - 581.5 kJ/kg ln k0 - 14.51 n (order) - 0.63 EA - 96.74 kJ/mol -159- Table 6.6 summarises some of the DSC data recorded for a number of different samples of skim milk powder. 1;b11_§;§ Calorimetric results for skim milk powders :9 AL _J.n_lso .11 (kJ mol) (kJ/kg) 96.74 581.5 14.51 0.63 * 83.4 421.6 11.97 0.6 85.53 482.48 12.42 0.58 87.52 530.61 12.55 0.61 81.69 445.73 11.62 0.6 ( * Avonmore product ) One clear point from this sample set of DSC data is that the self-heating reaction in the milk powder sample is not a zero order reaction. This is obvious in the context of the very small sample size used in the DSC sample cell. The above results would suggest that the reaction is nearer to first order in product concentration. As the results in Figure 4.6 show, there are two major peaks occurring in the thermal reaction of the powder. The range of operation allowed a complete kinetic evaluation of the first of these peaks. This was for the decomposition of the milk protein. The figures in Table 6.6 refer in the main to this reaction. The second reaction corresponds to the exothermic reaction involving lactose, beginning around 500'C. The evaluation of this peak was outside the range of the calorimeter. However, when the powder is at this temperature it is already in the ignition phase. Hence the data in Table 6.6 refer to -160- the exothermic reaction which instigates the powder fire. This is the origin of the fire trigger and hence the reaction data needed for the predictive studies is as presented in the table. The kinetic data needed for input into the F.E. model are the values for EA and Qf for the Avonmore skim milk powder. EA is 96.74 kJ/mol (see Table 6.6). The heat of combustion term,Q , is expressed per unit volume; hence the Exothermic Reaction energy term H must be converted from J kg“1 to J m"3 . The milk powder sample density is determined as previously outlined in Section 4.2.2. Thus, the Heat of Combustion of the sample is : Q = p AH - 600 x 581.50 x 103 a 3.49 x 108 J m'3 Including the frequency factor term , f (or kg in the DSC print-out), gives: Qf - 3.49 x 108 x 2.003 x 10 5 = 6.99 x 1014 w m'3 Table 6.7 allows a comparison of the DSC results with oven-derived values and with data published by Beever (1984 ) -161- 13213_§&1 Kinetic data for milk powder, derived by traditional and DSC techniques. ma 8.0.03.1: 12.8.9 E; (kJ/kg mol) 79.3 78.98 82-97 of (wm'3) 2.11 x 1013. 1.572 x 1013 7 x 1014. (* Avonmore ) The results indicate that the Activation Energy, EA, is given a higher value by the DSC than that predicted using the oven technique of Section 6.1. Beever's results were also obtained using the traditional oven techniques and are in line with the experimental values recorded here . In the particular case of the Avonmore skim milk, the EA value of 97 kJ/kg mol is significantly higher. This powder also shows as having a higher than average Heat of Combustion term. In terms of ignition then, while it takes a higher threshold temperature to initiate the self- heating reaction, it may subsequently proceed quite rapidly to a runaway reaction/ignition. Walker et al. (1983) also refer to this compensation effect whereby EA and the pre-exponential constant may be shown to be inter- related. The general effects of higher or lower EA values are considered further below when performing the sensitivity analysis. -162- £12.5139LAIIQE_B§§ELIE A FORTRAN program implements the finite element solution procedure outlined in Chapter 5. The solution procedure allows variation of time and space increments. No stability problems are encountered over the range of time/space increments used (e.g. time increments from 0.25 sec to 0.25 hours were used; space increments from 0.01 mm up) Small improvements in accuracy are theoretically possible with shorter time or distance steps. Methods of saving on computer time, such as outlined in the theoretical study of Anderson et a1. (1974), can not be used due to practicalities associated with data from 'real' product. In the idealised case used by Anderson et al., the time step was lengthened whenever the change in temperature was not significantly speeding up the simulation exercise. This procedure is not practical when dealing with real product and where the location of a center of ignition may also be important. Thus, while temperatures at or near surface nodes may have equillibrated, internal nodes may still be in the heating- up phase. A similar situation arises when self-heating commences, typically at or near the center of a sample. The onset of ignition, and hence the MIT, can only be determined by carefully monitoring the temperature across the sample using a fine grid and a small time increment. As the program is interactive, it is possible to institute -163- some savings in time by altering the selected time step. Thus, in zeroing in on MIT values , large time steps may be initially used to get a rough estimate of the location of the MIT. Finer selection of the MIT may use smaller time steps. The simulation program, listed in Appendix III, is designed to allow variation of the following parameters : Enzixgnmgntgl Ambient temperature Surface heat transfer coefficient fignplg Initial temperature Radius Density Specific heat Thermal conductivity Beastien.Bate.L.Einetie.nate Heat of reaction / combustion Arrhenius pre-exponential constant Activation energy .E1111_DIIBIL£D_§IHELAIIQE_BEH To illustrate the performance of the simulation model, a typical set of milk powder data is assembled. The results are shown for the powders as they undergo various stages of non-critical and critical self-heating. A typical skim milk powder sample has the following composition : -164- on e _& Moisture 4.25 Fat 1.49 Protein 33.46 Ash 6.80 Lactose 54.00 ( Source : Manufacturer's specification, Avonmore ) £111111_§2££IIIQ_EILI Equation [ 3.52] was used to determine the powder specific heat based on the compositional data. Thus, in the case of the powder sample, Equation [ 3.52 ] becomes: Cp = 1.424 X 0.54 + 1.549 X 0.3346 + 1.675 X 0.0149 + 0.837 X 0.068 + 4.187 X 0.0425 - 1.547 kJ/kg K i12111Z_AQIIILIIQE_IEI8QI_LED_IIBI_QI_£QHE!§119E The kinetic reaction rate data were determined from the plot based on the MIT results reported in Section 6.1 and listed as Table 6.3 : EA 8 79.316 kJ/kg mole Qf a 2.11 x 1013 W/m3 EA and Qf values are typical of the range encountered for milk powders. Similar values of 78.983 kJ / kg mole and -165- 1.572 x 10 13 W / m3, respectively, were obtained for the powder MITs reported by Beever (1984) (see Table 6.7). For the set of powders listed in Section 6.3.1, a Time / Temperature profile was developed using the simulation program. The results of the profile are shown in Figure 6.6 for six nodes, four equidistant internal nodes, the sphere surface and the center-point. MIT values may be pinpointed from plots such as those in Fig. 6.7. For the 4' sphere, a theoretical ignition is predicted at 136'C: at an ambient temperature of 135'C the powder temperature reaches a non-increasing equilibrium. The accuracy of the traditional oven method may only be to +/- 10'C ( Synnott et al., 1986). Predicting to l'C is quite accurate both in this context and when extrapolating results in an industrial environment. Clearly, the Finite Element model can be used to ascertain a particular MIT to whatever level of accuracy desired. but. better than 1'C is of little engineering value. The Finite Element program is established with the relevant. data on ‘the jpowder' characteristics ( 'thermal conductivity, specific heat, density) and the kinetics of -166- Temperature ('C) r 1500 I 1000 Time ( mine ) Figure 6.6 Simulated Temperature Profile for 4" Sphere -167- Temp. ('C) 30°F ml- 400- 136'C m. 135'C P o A l A L A l 0 500 1M 1500 Time (mine) Fig. 6.7 Predicting MIT Using Oven Date. 1678 the self-ignition reaction ( Activation Energy, Heat of Combustion.x.Frequency Factor). The simulation process is initiated by entering data on the sample size, ambient temperature, iteration step size and the number of elements required in the Finite Element grid. The time / temperature information is recorded for each node after each time interval . Subsequently, the ambient temperature is varied to find the lowest temperature at which the powder ignites, within an accuracy of 1'C. Thus Fig. 6.7 shows non-ignition of a sample at 135'C. When the ambient temperature is raised to 136'C, ignition takes place, as the exothermic reaction takes over. One of the advantages of the simulation process is that long runs can readily be simulated. This ensures that situations involving a slow, smouldering fire (which may, in typical industrial storage conditions allow a fire a lengthy 'incubation period') , can be simulated without recourse to long oven tests. Besides, the experimental tests may prove to be inconclusive, due to practical limitations on powder/reactant mass etc. W The aim of the simulation program is to simulate the process of self-heating, leading potentially to self- ~168- ignition, in the milk powders under investigation. The validity of the model is illustrated in this section. The model was first used to establish the Minimum Ignition Temperature values for the samples used in the experimental work. The values for the Activation Energy, the heat of combustion, etc. were taken from Table 6.5. Figure 6.7 shows the MIT value for the nominal 4" sphere to be 136°C, with a non-ignition at 135°C. The plot may be compared with Figure 4.5, the experimental plot for the 4" sphere of Avonmore skim milk powder. Figure 4.5 shows where an oven temperature of 137.27'C results in non- ignition, with an ignition at 138.9'C, Table 6.3. The experimental MIT is 138.09'C, compared to the simulated predicted temperature of 136°C shown in Fig. 6.7. Table 6.8 compares the predicted and experimental MIT values for each of the standard sphere sizes. Iihll_£1§ Oven a Model MIT values compared. Heminal Experimental Model Lesser Diameter 811 MIT Ignition (Inches) (°C) (°C) (°C) 1.5 172.64 168 171.46 2 161.18 159 162 3 143.91 144 146 4 138.08 136 138.9 The Experimental_fi11 is the standard 'average' MIT value, the ngg§t_1gni;ign is the lowest experimental temperature at which ignition was recorded. '169- The simulation model is within 1°C of the oven MIT for the 3" sphere, and is within 2°C of the oven value for the 2" and 4" samples. To the nearest °C, the variation of the predicted value ranges from 0°C for the 3" sample diameter to 5°C for the 1.5" sphere. As Table 6.8 shows, the model MIT is generally 2°C lower than that determined by the oven technique. Thus, the model, based on a single Arrhenius reaction term, and a knowledge of the product composition and the environmental conditions, gives a conservative estimate of the MIT. From a safety point of view this is fortunate. Some interesting results arise when the lowest oven ignition temperature is used as the simulated oven ambient temperature. Thus, for the 1.5'° sphere the lowest ignition temperature was found at an oven temperature of 173.82'C. Using this as the ambient temperature in the model, produces the time/temperature profile shown in Fig. 6.8. While the ignition-induction time is longer for the model, the peak temperature for the oven test is reached almost simultaneously by the simulated profile. Time to Ignition in these plots is c.140 minutes for both the simulated and experimental cases. Fig. 6.9 for the 4" sphere shows the same phenomenon for the time/temperature profile at 138.9'C, (i.e. at the 'lowest' ignition temperature for this sphere) The profiles of Figs. 6.8 and 6.9 have strong ~170- 1000 - m u- .3 a 500 ’ e I I— - . 400 - Oven a Model .v" x __.p"" 2m 14’5””! 7'" . ,4” .'e’.;" 0 A l . J 0 100 200 Time (mine) Fig. 6.8 Oven 8. Model Profile at Lowest Oven Ignition Temperature (1.5" Sphere). '171- implications for the model, for its usefulness and validity. The simulated profile does not exactly follow the oven profile in its entirety. But, with respect to the main accelerating self-heating reaction, the model is accurate. This is further borne out when a similar comparison is made with the 2" and 3" test spheres, as shown in Figs. 6.10 and 6.11 respectively. Here again the induction times for both experimental and simulated plots correlate well. Thus the model can safely and accurately be used to predict both MIT and Time to Ignition for the milk powder. Given the complex nature of milk powder, a number . of components are oxidising and producing heat as well as that described by Equation [3.2] used in this analysis, i.e., G = Q f e'EA/RTamb [3.2] As the sample graph in Figure 4.6 shows, DSC analysis identifies two main exothermic reactions taking place in the milk powder. The first reaction is due to the decomposition of protein which begins about 220°C, followed by a much larger exothermic reaction above 400°C caused by the decomposition of lactose. As the compositional breakdown shows , protein and lactose are -172- Oven Temp. (°C) Model 0 200 400 600 800 1000 Tlme (mine) Fig. 6.9 Oven a Model Profile at Lowest Oven Ignition Temperature (4" Sphere). 1728 m .- 1; Oven 60° ’ < ‘. 53. i \I. i I. d 400 " f . I § . . Model b I ' ‘- i- f. :41 m I- 0 - . . 0 200 400 600 800 1000 Time (mins) Fig. 6.10 Oven and Model Profile Oven Ignition Temperature (3"). at Lowest -173- Temp ('C) 1(00" m I. 1"} i f E ‘TVEN ‘ 5. m .- I . 400 - g MODEL 400 6m 8m Tlme (mine) Fig. 6.11 Oven 81 Model Profile at Lowest Oven Ignition Temperature (2"). 173a the two principal components in the skim milk powder. The close agreement between the experimental temperature history at the lowest ignition temperature and the predicted profile indicates that the model is well in line with 'the jprime self-heating reaction. It is therefore concluded that it can be used with confidence to predict and study the phenomenon. When using it to predict the MIT values, the model errs, at least with respect to the oven values, on the side of safety. A more detailed chemical / kinetic analysis is required to isolate the reactions occurring at the higher temperatures, where the oven profile temperatures increase faster due to secondary reactions. Thus at 138.9°C, the simulated temperature rise occurs at c. 700 mins. The sample in the oven shows an increasing temperature (although at a slow rate) after c. 450 mins. The faster induction time , not forecast by the model, is compensated by the lower simulated. MIT values. Thus, the general trend is for the ‘model to predict ignition at a lower MIT value than. is actually observed in the oven tests. A number of factors, experimental and theoretical contribute to the divergence. As the model is currently programmed, a fixed value of thermal conductivity is used throughout the simulation. The thermal conductivity of the powder actually increases with -174- increased temperature. Experimental k-values reported by Martin (1987) on skim milk powder are of the order of 0.1 W/mK at 100°C. Thermal conductivity values at high temperature are in fact difficult to determine. Martin's results were obtained using a guarded-plate technique. The increase in thermal conductivity as the powder heats up to ambient and above, means that the heat may be transferred faster from the hot sphere center to the relatively cooler surface area. This heat transfer allows the powder to tolerate the self-heating for a longer time, thereby postponing ignition. Conversely, a higher ambient temperature is needed to initiate eventual ignition, i.e. the powder has a higher experimental MIT than the predicted MIT. This probably accounts for some of the model / oven divergence in Table 6.7. A Certain experimental difficulties were also encountered in trying to get more accurate control of the oven temperature. While the temperature controller maintained a steady temperature within the oven, the correlation between the controller setting and the ambient temperature within the oven did not permit fine resolution of the temperature set point. é11_DflQ_A_QEEE_£_§IH!LAIIQE_QQHILBI§QE As seen from the results shown in Table 6.7 , the DSC -175- values for the Activation Energy and the pre-exponential and Heat of Combustion term, Qf, are significantly different from the oven-derived values. For the Avonmore skim milk powder, the 4" sphere has a DSC EA value of 96.74 kJ /kg mol compared to an 'oven-derived' value of 79.3 kJ / kg mol. The implications of this are clear: the higher the Activation Energy, the more external heat is required to initiate the runaway thermal reaction. In effect, the high EA value leads to a higher MIT or longer Time to Ignition. This effect is offset in some degree by the higher value of Qf (i.e. 7 x 1014 vs. 2.11 X 1013 W m'3). Thus, when the reaction is in progress, the rate of heat generation is greatly increased. Figure 6.12 illustrates how the Finite Element model may be used to predict a MIT value for the 4" sphere using the kinetic data derived from the DSC results. Simulating the self-ignition phenomenon with the higher value of EA and the higher value of Qf, leads to a MIT of 159°C, in comparison to the values of 136°C and 138°C obtained using the oven derived kinetic data and the direct oven MIT value, respectively. Table 6.9 extends this comparison to the other standard sphere sizes. -176- ('C) 8 see as Temp. Fig , 159.0 g P /. 158.0 100 -, J I 0 . 1 L 1 . 1 0 1000 2000 3000 Time (mins) . 6.12 Predicting MIT Using DSC Data for a 'Standard' 4" Sphere 0f Skim Milk Powder. -177- 13h1g_§;2 Comparison of MIT(°C) values for different size spheres using different techniques. we MILIJIEJ. 12 32 22 1.5." OVEN 138 144 161 173 MODEL/OVEN 136 144 159 168 MODEL/DSC 159 167 179 188 While the oven and model/oven values correlate well, the MODEL/DSC MIT values are significantly different. This is due to the threshold action of the high EA value measured by the DSC. The DSC also predicts higher levels of heat generation. As the model and oven profiles show, once the reaction's energy threshold has been breached the runaway reaction. takes over, and. the reaction ‘proceeds almost instantaneously to ignition (Figure 6.12). Since predicting this ignition point is the main thrust of this study, employment of the DSC is warrants further study since it appears to assume an inaccurate representation of the reaction (up to ignition). £1511_DEQ_LQQ!BLQI The DSC printout in Figure 4.6 shows the 90% confidence limits on which the DSC calculations / regression are based. This leaves a potential error of +/-10%. Applying this error criterion to the DSC MIT value in Table 6.9, gives a potential error of approximately. 16°C (i.e. a -178- 'corrected' potential OVEN/DSC MIT of 143°C, only 5°C higher than the MIT value determined by the conventional oven test). Viewed in this light, the DSC results are a further verification of the quality of the model. This also favours future use of the model in conjunction with Differential Scanning Calorimetry, and may lead ultimately to the replacement of the time-consuming traditional convection-oven or hot-plate techniques. £1§IZ_QQEQL!§IQE§ As the DSC data shows, the DSC analysis assumes a non-zero order model, while the F.E. model simulates a zero order system. Forcing non-zero order data into the model does account for the inherent inaccuracy in the DSC/MODEL results. This is an important factor in trying to correlate the three sets of results. Values of the Activation Energy and the Heat of Combustion X the Frequency Factor also make a significant contribution to the deviation between the OVEN and DSC results. Differential Scanning Calorimetry was used to further improve the overall simulation/prediction procedure. When determining the highest non-ignition with a convection oven, several long test times may be necessary, of up to -179- 2000 minutes in some cases (IDF, 1987). In contrast, a complete set of product kinetic data may be obtained from the DSC in less than three hours. Employing the DSC data in the model introduces some experimental/operational difficulties. In particular, there are difficulties in terms of interpretation of the results. One significant factor is the question of the order of the reaction. The model assumes a zero-order reaction - an assumption safely made for the industrial situation; it is of doubtful merit for use in the convection oven tests and is not true for the DSC sample cell. The reaction is much closer to to a first order than a zero order case. The non-zero order of reaction term can only be used in the model if the experimental record includes time/temperature and time/mass data . This adds a further complexity to both the model and the experimental verification. It also renders the model more difficult to use outside a research laboratory. When instrument and regression accuracies are accounted for, ignoring the non- zero order factor does not lead to significant errors. Hence, the combination of the model and the DSC has potential for rapid analysis of powder combustibility. Further work is needed to establish the error bounds and applicability for the DSC kinetic data. -180- §1£_££E§III!IIX_LEBLX§I§ After establishing the quality of the model with a heat generation term, the model can be employed to conduct a sensitivity analysis on the factors influencing the combustion / ignition of milk powders. Table 6.10 shows the range of parameters used in the various simulation 111118. 25b11_§;10 Simulation parameter range. Bannister 5mm Bangs Activation Energy EA 50 - 100 kJ/mol Sphere Radius r0 1 - 12 cm Ambient Temperature Tamb 150 - 400 °C Surface Heat Transfer Coefficient h 5 - 25 W/m2K The range of values was selected to center on a typical set of milk powder parameters shown in Table 6.11. The temperatures selected are typical of powder manufacturing , conveying or storage conditions (of. Figure 1.1) and are within the temperature range which could be achieved using the test oven, Section 4.2. The 'standard' values for the simulation runs are listed in Table 6.11. -181- 13p1g_§.11 'Standard' Simulation Data Activation Energy Sphere radius Ambient Temperature H.T. Coefficient Density 79 kJ / mol 5.08 cm 195 'c 16.7 w / mzK 600 kg / m3 Time / Temperature profiles were obtained for the 'standard' conditions listed in Table 6.11 for different sphere radii. Minimum Ignition Temperatures were obtained, to an accuracy of 1°C, for each of the radii. Thus for the 4 cm radius, ignition was found to occur at 153°C. At 152°C, though the sample temperature increased [ to 155°C at the surface and 178°C at the sphere center, after over 4000 mins ], there was no ignition / runaway reaction. In the case of non-ignition, this represents a state of equillibrium between the heat generation and the surface heat loss ( cf. Fig. 1.2 ). Table 6.12 shows the simulated MIT value as a function of the sphere radius. 13213_§;1z MIT vs Sphere Radius. Basins :11 (cm) (°C) 1 207 2 180 4 153 5.08 145 6 138 8 129 12 117 -182- The data in Table 6.12 is plotted in Figure 6.13. The results clearly show an exponential increase in MIT with decreasing sphere radius. For large spheres , ignition occurs at relatively low ambient temperatures, e.g. for a 12 cm radius, the predicted MIT is 117°C. As the sphere heats up in the convective air-stream, the larger the sphere the longer it takes for the heat to reach the surface to be dissipated. This causes a build-up of heat internally, leading to a further increase in the rate of heat generation: this in turn causes the temperature to increase further, leading ultimately to thermal runaway and ignition. Smaller spheres need much higher temperatures to ignite as they can more easily dissipate the heat generated internally. Thus, the trend shown in Figure 6.13 is in agreement with normal heat transfer considerations. The influence of sample size on MIT also agrees with experimental observations as summarised in Table 6.3. The oven results qouted show MIT decreasing from 173°C to 138°C as the sample sphere radius increases from 1.5" through to 4". The IDF summary document (IDF,1987) confirms this trend. The IDF proposes as a general guideline a 15°C decrease in MIT for the doubling of sample radius for the general category of dairy powders. The results of Table 6.12 would suggest a decrease of 20°C '183- 200 180 ('C) 160 M.|.T. 140 120 100 Figure 6.13 Simulated MIT vs. 2 4 6 8 10 12 14 Radius (cm) Radius for Milk Powder. -184- or higher in the particular case of skim milk powder. The trend shown in Figure 6.10 is also in agreement with the work of Boddington et a1. (1981) who predicted that under the Semenov conditions of low Biot number ignition is less likely as reactant mass size is reduced. Boddington et al's case of negligible internal thermal resistance is a subset of the solutions presented by the simulation process here. £151112_IIHI_IQ_I§!IIIQE While the MIT decreases with increased critical dimension/ sphere radius, increasing the radius of the sample has the opposite effect on the Time to Ignition (TtI). To investigate this, the TtI was predicted for a range of sphere sizes, .with a fixed ambient temperature of 195°C keeping the other 'standard' parameters as listed in Table 6.11. The results of this analysis are shown in Table 6.13 and plotted in Figure 6.14. -185- Time to Ignition (mine) Fig 300 I- 200 . 100 - o A 1 1 1 A l A 1 A L L l A l 0 2 4 6 8 10 12 14 Radius (cm) 6.14 Simulated Time to Ignition vs Radius. -186- ng1g_§;1; TtI vs. Sphere Radius .131 (mins) 73 93 122 156 180 202 224 256 283 WEE HHmQChUlb I O 00 NO From Figure 6.14 it evident that the Time to Ignition increases steadily as the radius increases. The relationship is approximately linear over the mid-range of radii , with the rate of increase tailing off slightly at greater radii. This happens since with increasing critical dimension it takes longer for the heat to transfer through to the center of the sphere. Thus, while the MIT is lower for larger spheres, the 'induction' time is significantly longer. This is borne out by experimental results which show that the large radius samples have longer heating-up times. The combination of low MIT and high TtI poses serious implications in commercial installations as relatively low ambient / storage temperatures may, over a protracted period of time , create critical conditions which leads to self-ignition of powder deposits or powder stores, (Avonmore, 1987). -187- As the results show, no ignition was recorded for the 1 cm sphere as the MIT for this sphere ( 207°C, Table 6.12 ), is above the ambient temperature of 195°C employed in this analysis. 41'!” " .‘1‘ TI. '1 diff." '1 91°11‘08 '1 , 0 Comparing the three sets of results obtained and summarised in Table 6.9, the importance of Activation Energy has been referred to already ( see also Table 6.7). The Activation Energy is varied in the simulation exercise to assess the influence on the combustibility of milk powder. The of value is held constant. In practice, following the sample calculation outlined in Section 6.3, an increased value of EA as obtained from the plot of ln ( Schch/roz) vs. l/Tcr, Figure 3.1, leads to a decreased value of Qf. Thus, a higher EA derived from such a plot slows down the beginning of the exothermic reaction: it also signals less heat generation at the outset of the process and hence a further 'postponement' of ignition. In the data recorded, Qf is held constant (i.e. Qf is independent of a variation in EA, the parameter under study). The results are recorded in Table 6.14, and are shown in graphical form in Figure 6.15. -188- Fig. 10X)- 5m)- Tlme to Ignition (mine) r 40 50 60 70 80 Activation Energy (kJ/mel) 6.15 Simulated Time to Ignition -189- V I 90 vs Activation Energy. 6.1 TtI vs. Activation Energy. W .211 (kJ / mol) (mins) 50 1.8 60 9 70 41 75 91 79 156 85 319 88 1256 A number of salient points emerge from the plot of TtI vs. Activation Energy. Firstly, it is clear that an Activation Energy threshold exists above which ignition does not take place. Thus, if the EA value for a particular exothermic reaction is too high for the energy content of the product / environment, the reaction will not proceed. In the case of the milk powder in Table 6.11, placed in the 16.7 W / m2 K air stream at 195°C, the threshold is approximately 90 kJ / mol. No ignition occurs at or above this value of EA. At the other end of the scale, when the EA is below a certain level, the reaction is practically instantaneous. Thus, for the reference conditions,at an Activation Energy of 50 kJ/mol the reaction accelerates in less than two minutes. Thus, it is very important to establish the value of the kinetic parameter when evaluating the level of risk in manufacturing or handling milk powders. As the experimental results already referred to indicate, a high. Activation Energy is often associated. with an -190- increased value for Heat of Combustion. This is shown in Tables 6.6 and 6.7. This compensatory effect was also proposed in the theoretical calorimeter model of Walker et al. (1983) Simulated results were obtained for the effect of ambient test temperature on the ignition parameters. A range of temperature is tested to establish the time required to ignite the 'standard' sample sphere. In the event of non- ignition, the analysis indicates the threshoLd condition for ignition. Table 6.15 and Figure 6.16, respectively, present the numerical and graphical results of this investigation. IAEL§_§&1§ Simulated Ambient Temperature vs TtI. Annient.xeneerature .111 (°C) (mins) 145 1225 150 530 175 219 195 156 225 90 250 58 300 29 350 17 400 11 -191- 13X)- 100:- 5a)“ Tlme to ignition (mine) 100 200 300 400 500 Ambient Temperature ('C) 6.16 Simulated Time to Ignition vs Ambient Temperature. -192- As with the plot of EA vs. TtI, a variation of TtI with ambient temperature has two distinct elements. No ignition is possible below the MIT temperature: as the ambient falls towards this temperature, the time to ignition begins to increase. As illustrated in Figure 4.5, the MIT for the standard sample is 138°C. From Table 6.15 it can be seen that at 150°C it requires almost 530 minutes to ignite the sample. At the other end of the scale, ignition is fast, requiring only 17 mins at 350°C, and 11 mins at 400°C. The exponential reduction in TtI with increasing temperature shows that above 300°C case, ignition is rapid. Thus, in an industrial environment, situation, short-term exposure to a high temperature source (e.g. a welding torch flame being applied to the external surface of a pipe or a silo) is enough to initiate a critical ignition situation. Reference to the experimental data collected confirms this trend. Thus, a difference of just 4°C in oven ambient temperature makes a difference of 100 mins induction time for the 4" sphere whose experimental temperature profiles are shown in Figure 4.6. As already noted, these relatively long induction times are associated with temperatures near the MIT. As temperatures increase significantly beyond MIT, induction times tend to merge and for all such temperatures as the exothermic reaction takes off very rapidly. -193- The effective surface heat transfer coefficient, h, is varied over the range of values shown in Table 6.10 to monitor the influence on Time to Ignition by h. Industrially, this value will vary depending on the location of the powder in the manufacturing or conveying system, Figure 1.1. A sample being heated up internally, i.e. undergoing an exothermic reaction, will be dependant on the ambient airstream to effectively cool the surface and thus dissipate the generated heat. Thus the effective h-value plays a dual role in the self-heating process. Initially it acts as the heat source supplying the thermal energy to initiate the self-heating exothermic reaction(s) within the powder. Thereafter, as the internal temperature of the sample rises above ambient, non-ignition depends on the ability of the ambient airstream to function as a heat sink for the generated heat. The surface heat transfer characteristics, combined with the size of the temperature differential dictate the success of this heat dissipation. The role of the h-value is analysed in this part of the sensitivity analysis. The simulated results are tabulated in Table 6.16 and shown graphically in Figure 6.17. -194- 220 _ g 200 ~ E v 4’ C 2 E 180 ~ 2 3 ,, I E 160 - p. 140 5 l . 4 5 ' 0 10 20 30 Surface Heat Transfer Coefficient (Wlm K) Fig. 6.17 Simulated Time to Ignition vs 'h' - Value. -195- 13§L§_§.;§ Heat Transfer Coefficient vs TtI 94.991251211141111 .1311 (W/m K) (mins) 5 210 7.5 180 10 168 16.7 156 20 153 25 149 Given the nature of the surface heat transfer mechanism, the results in Figure 6.17 are as expected. For low values of the heat transfer coefficient, the sphere surface is slow to absorb heat from the environment, and hence self- ignition is delayed. When 'h' is too low, the surface thermal resistance prevents ignition as the sample temperature does not heat up fast enough or high enough to initiate a runaway exothermic reaction. By contrast, when 'h' exceeds approximately 15 W / m2 K, in Figure 6.17, the surface heat transfer has reached its optimum level. In these cases heat transfer is quite rapid and the sample reaches ignition temperatures very quickly. While the increased surface heat transfer aids in dissipating excess heat from the surface, which during the course of the exothermic reaction is at a higher temperature than ambient, the finite nature of the internal thermal resistance mitigates against this. Thus the high surface heat transfer coefficient effectively initiates ignition once Activation Energy and Ambient Temperature are in line with the conditions outlined in Sections 6.4.2 and 6.4.3. -196- Extrapolating this curve to an infinite 'h' value, the original Frank-Kamenetskii assumption, yields a TtI in this instance of approximately 145 mins. This allows the model to readily simulate the typical heat transfer conditions allowed for in the literature. The plot in Figure 6.17 emphasises the error inherent in the Frank- Kamenetskii heat transfer assumption, particularly in relation to situations of low surface heat transfer coefficient e.g. storage conditions. Strict application of this assumption would result in a significantly shorter TtI prediction and also a lower MIT value. -197- -198- 11.52HHABX This study set out to develop a simulation process to model the self-heating/self-ignition phenomena which are a major source of concern and of significant financial loss, to milk powder manufacturers. Firstly, the model was shown to be a true representation of the phenomena. Next it was used to identify the chief environmental factors which affect self-ignition. On the experimental side, the use of the Differential Scanning Calorimetry technique to replace the conventional oven tests for supplying raw kinetic data to the model was studied. The theory of self-ignition has been firmly rooted in the Frank-Kamenetskii approach for more than half a century. Successive authors have contented themselves with peripheral, empirical embellishments. of ‘the Frank- Kamenetskii model. The F-K model has as its central tenet the assumption that there is negligible surface resistance to heat transfer, and that the reactant mass is limitless. The present work casts doubt on both assumptions. Finite surface thermal resistance may pertain to many potentially explosive industrial storage/handling conditions. A variable surface 'h' value affects the product Minimum Ignition Temperature and the Time to Ignition. With regard to the reactant concentration, laboratory tests , in either convection or DSC oven, need to be interpreted in -199- the context of scale-up. The DSC with its highly restricted sample size cannot replicate zero-order conditions. The simulation model presents a stable and versatile means of predicting the likelihood of spontaneous ignition during milk powder manufacturing. The model is a major step forward as it makes no assumptions on the crucial question of either internal or surface resistance to heat transfer. Two types of data are needed to employ the model : [1] Environmental and physical data, in particular product sample size and the surface heat transfer conditions. [2] Product property data - thermo/physical data and kinetic data. Data of this first type may be known or can be obtained by basic site investigations and/or laboratory-scale modelling. Product property data (i.e. kinetic data) is not available presently in the literaturefor the majority of the broad range of commercial milk powders. Thus, it must be determined independently. Two forms of determination may be used: -200- [1] Traditional oven techniques, and [2] Differential Scanning Calorimetry. The oven-based kinetic data determination gives the better correlation between the model and the oven MIT and the general time/temperature profiles. Thus, after establishing a powder's kinetic data via a number of oven tests at various temperatures, the other temperatures, sample sizes, andheat transfer conditions may be simulated with confidence. This represents a major step forward in determining combustion parameters for milk powders. Results show that [1]Product Activation Energy and Ambient Temperature together dictate a product's liability to ignite. [2]Effective powder handling can prevent self-ignition. When large clumps of powder form in a dryer or silo the powder MIT is lowered and the Time to Ignition is increased. Thus, powder in storage may ignite, even at the relatively lower storage tempereatures, after a lengthy period of time. The use of a conditioned air powder conveying system will reduce this risk considerably. -201- E12_EEQEBfiIIQE§_EQB_ZEIEBI_fiI!DI 1. Conduct a more detailed experimental and theoretical investigation on the effect of reactant consumption on critical combustion parameters. As the DSC results show, the exothermic reaction central to the model is not a zero-order reaction. This is true where, e.g., in either the convection oven or the DSC furnace, small samples are studied. To properly establish the scale-up factors, the simulation should employ the general heat generation term G - Q f Ci" e‘EA/RT 2. Develop specific applications of the model to simulate known areas of high risk in powder plants, e.g. dryer 'hot spots', storage siloes, fluidised beds. Specific initial and boundary conditions should be used in the model to account for high-risk industrial situations. The simulation results should be correlated with recorded data on known powder fires to verify the applicability of the method. This would also greater acceptability for the technique within the dairy industry. 3. Select the kinetic and thermal property data employed in the model to simulate sensitive dryer 'start-up', conditions e.g. contaminated, overheated or rewetted -202- powder deposits in dryer. These have been implicated as hazardous operating conditions in a number of industrial fires, Synnott et al. (1986) and Beever (1984). This work would entail significant work in preparing specific batches of spray-dried powder, followed by detailed property determination. 4. Investiagte the effect of milk powder compositional factors, e.g. moisture, fat content on self-heating rates. Isolation of the milk components causing the significant exothermic reactions would facilitate effective classification of powders in terms of ignitability. High temperature calorimetry would assist in this process. 5. Extend the program applications to non-dairy and, possibly, non-food products. While the samples used here were all dairy powders, the simulation is based on first principles and so will work with similar effect independent of the product type. To date models used for example by Walker et al. (1983) for wool have been empirical and were limited to the extreme boundary conditions of either the Frank-Kamenetskii or Semenov models. This would be an important contribution to studies in fire prevention. 6. Develop a model application to simulate start-up losses in spray dryers, e.g. the production of 'first run' burnt '203- particles. Loss prevention and minimisation are becoming of major interest to powder manufacturers. A finite time is needed to correctly balance the flows of heat and product into the dryer. In spite of much research, dryer start-up and control are not yet automated. For bacteriological as well as fire safety reasons, dryers are started up on the hotside, i.e. with a view to avoiding the transport of partially dried or damp product into the powder conveying and storage areas. This type of powder can lead to contamination, and to the formation of moist, self-heating clumps of powder. The result of the 'overheating' of the first batch of product through the system is the presence of burnt particles in the system which must be retrieved and discarded. The present model may be used to simulate part of the heat balance during this critical start-up phase, and hence minimises the production of the unwanted byproducts. 7. Source or otherwise determine the high temperature thermal properties of powders. Their influence on the combustion parameters may then be evaluated. The results indicate that increasing the thermal conductivity leads to a higher effective MIT. Thus, a built-in calculation / database of conductivity vstemperature will improve the model predictions by one or two degrees Celsius. -204- 8. Determine the effect on the MIT of the initial temperature of the sample. This is of particular relevance where the powder is precooled between the drier exit and the entrance into the storage siloes. This is a common practice in the dairy industry as it guards against potentially hazardous temperatures on storage. In some climatic conditions (i.e. where humidity is high) the conveying air is precooled to ensure the dry powder is not put in contact with moist air prior to storage. -205- APPENDICES APPENDIX 1 Derivation of Element Matrices To establish the element matrices, matrix algebra manipulation is used. A number of definitions are needed. The temperature 'within. an element. is assumed. to vary linearly' with distance from. the nodes: i.e., for the general "e" element : T(e) = ci(e) + cj (e) A-l where ci(e) and Cj(e) may be determined for a particular element from the known nodal temperatures Ti and Tj. Thus: Ti = ci(e) + Cj‘e) ri and Tj = ci(e) + Cj(e) r5 , A-2 giving ci(e) = rj Ti - r1 Tj rj - r1 and A-3 '206- cj(e) g Tj - Ti rj - ri To facilitate the determination of the various element matrices, Equation A-l can be written in matrix form as : T(e) = [ 1 r ] [ci] (e) A-4 Ci Following the method of Myers (1971), the following matrix definitions are useful : pT = [ 1 r ] A-S c = Ci (e) Cj A-6 Thus Equation A-4 can be stated as: T(e) a pT C(e) A-7 Equations A-2 , can also be written in matrix form : T1 = 1 ri ci (e) Tj l rj Cj A-8 -207- The matrix on the left side of Equation A-8 is 1(9), the matrix of nodal temperatures . Equation A-8 yields a new matrix, defined as : P(e) = 1 r1 1 rj A-9 Thus, 1(9) = P(e) C dV 2 4 w r2 dr -209- => IRIS) = 2 «I rj kr(e) r2 (ell—3.0; ri (dr) => = 2 r kr(e)J rj r2 (d1) 2 dr ri (dr) A-19 kr(e) ( written simply as k hereafter) is taken outside the integral sign as it is assumed to be constant over the element. This is a reasonable assumption, provided the element is small enough that the temperature / thermal conductivity does not vary within the element. Using matrix notation, Ik(e) may be expressed as : Ik(e) = 2 I k I rj r2[ pTrR(e) I ]2 dR ri A-ZO Differentiating d1 (9) = 4 r r' T R (e) d T RT r2 d dike) JriJ(P r I ) (P(£) ) r d1 A-Zl Let pTr R be defined as aT - [ a1 a2 ] (constants independent of 1(9) and r ) Then : aT = To I [ a1 a2 1 1 T3 = alTi + asz A-22 0(aTI) a1 d1 = = (aT)T = a a2 A-23 -210- géeirsri = (9T: RIT fl: d1K r1 = 4 I = 4 1 = 4 r k RT r1 = 4 r k RT rj ri 0 = 4 r k RT 0 Premultiplying by RT 161 rj -1 d2 rij -r1 1 0 -1 r' O l -1 1 -1 = X -1 1 b .J where X= 4 t k (rj3 - r13) 3 rij3 and rij - rj - ri Therefore K(e) . 4 I k (rj3 - r13 ) 3 rijz -211- l 0 1 0 1 the element conduction matrix is A-24 4 x k rIIpTrRTI(pTrRITr2dr k l r3(PTrR)T(PTrRI)r2dr ri k I rj RT pr pTr RI r2 dr ri r1 prpTr r2 dr R I J [0 1]r2drRI (rj3 -ri3) R T 3 J as indicated 0 0 11111315111212 0 1 I 1 A1A1_QABLQIILEQI_EL§H£EI_HLIBIZ From Equation [5.6], the heat capacitance component in the 'functional' is Ic(e) = p c T (OT/0t) dV v v = (4/3) x r3,=> dV = 4 x r2 dr Ic(e) = 2 r I rj pc (2112 r2 dr ri (at) 2 1 pc g_ rj T2 rzdr at ri 2 r pcg_ rj (pT R(e)1(e))2 r2 dr at ri using the identity established by Equation A-14. Minimising with respect to temperature 1(9) gives : dIc(e) =4 1pcg_ rj(pTRI)(pTR)T r2dr d1 at ri A-ZG Equation A-26 follows from the differential rule set out in Equation A-23, following from the fact that both pT and R(e) are independent of temperature ( Equations A-5 and A- 9 ). Matrix algebra allows the integral term to be written in a more convenient form as : dIc - 4 1 p rj RTppTRI; r2 dr d1 ri (1; denotes differentiation w.r.t time, t ) = 4 1 pc RTI rj ppTr2 dr R I; ri ( R, Rt and I; are all independent of r and hence may be taken outside the integral sign ) [ rj r2 r3 = 4 1 pc RT dr R I; Jri r3 r4 J rj -1 rj3-ri3 rj4-ri4 = YY 3 4 R I; -r1 1 rj4-ri4 er-ri5 5 AI 4 5 L d ~212- ( Y! = 4 1 pc is defined to simplify the notation ) When the matrix multiplication is completed as shown, a 2 X 2 set of coefficients results defining' the general capacitance element matrix : 4 7 pc C11 C12 [cp(e)] = A-27 60 (ri-rj)2 621 622 where : C11 = 2rj5-20rj2ri3+30rjri4-12ri5 c1; = 3rj5-5rj4ri+5rjri4-3r15 c21 ‘ C12 c2; = 12rj5-30rj4ri+20rj3riZ-2r15 A-28 LA1I1_QQEEZQIIQE_£L§HEEI_HAIBIZ Matrix algebra may similarly be used to develop the convection element matrix. However, since convection only affects the final or Eth element, a simpler technique is employed here. From Equation [5.6] : 1(8) = 1/2 I h(T — Tamb)2 dS A-29 S Since h only has a non-zero value at the surface or Eth element, 1hE+1 = 1/2 J8h(TE+1 - Tamb)2 ds A-3O -213- Assuming that h is uniform ( i.e. constant ) across the surface and TE+1 is constant, the equation simplifies to . IhE+1 = 1/2 h (T311 - Tamb)2 IS dS A-31 But I dS = sphere surface area = 4 r R2 3 => IhE+1 = 2 x h R2(TE+1 - Tamb)2 A-32 Minimising dInEil = 4 r h R2(TE+1 - Tamb) d1 = 4 w h R2 TE+1 - 4 x h RZTamb 0 0 TE -47rhTamb u a a D’ 0 32 TE+1 A-33 ( written in matrix form ) -214- APPENDIX 2 DSC Programing Procedures 31211_§£B£§EIE§_H§IEQD The screening method is used for the initial thermoanalytical study of a substance. It is generally applied over a wide temperature range at a medium to high heating rate ( 10-50 K/min). Upon insertion of the sample, the following instructions are entered via the keyboard: INPUT SCREEN DISPLAY [ SCREEN ] Scan Parameters [ RUN ] Start Temp °C 50 : Input of Start Temp required [ RUN ] Rate K/min 10 : Heating rate required [ RUN ] End Temp °C 250 [ RUN ] Time Iso Min 0 [ RUN ] Plot cm 20 [ RUN ] Offset 80% [ RUN ] Pan type 1/2 1 [ RUN ] Limit mW 0 Evaluation parameters Screen [ RUN ] Dyn/Iso 1/2 0 [ RUN ] DSC °C 35 [ RUN ] Ident No. [ RUN ] Weight mg : Insert weight of the sample [ RUN ] Insert °C [ RUN ] Insert sample The DSC proceeds to increase the temperature at the specified rate. With the printer on-line, the DSC curve is printed out for the process. -215- 81212.5!ALEAILQE The DSC curve can be evaluated by the Processor using methods such as Integration or Kinetic . The evaluation component of the method is called up using the appropriate keyboard sequence. In addition, the Screening method has an evaluation function which allows a portion of the selected section of the curve to be plotted. Thus : DSC Screen 'C 45 Evaluate ] Evaluate I [ Screen ] Screen [ RUN ] Dyn/Iso 1/2 0 :Input 1 [ RUN ] Start 140 :Input required start Temperature [ RUN ] End 200 :Input required end Temperature [ RUN ] Baseline type 1 [ RUN 1 Plot cm 10 [ RUN ] Plot mode 1 [ RUN ] Weight mg [ RUN ] *** Calculating *** DSC Screen °C 50 A plot of the evaluated portion of the graph is then printed out. 3.2.1.Infaszafien Using the integration function, the heat flow of a dynamic or isothermal experiment is integrated allowing a direct calculation of the heat/enthalpy change in the observed reaction. The Integration input parameters are as listed. DSC MS 'C 35 [INTEG] Scan parameters [ RUN ] Start Temp °C 100 [ RUN ] Rate K/min 5 [ RUN ] End Temp °C 145 '216- Hf—IHHHf—I RUN RUN RUN RUN RUN RUN Evaluation RUN RUN RUN RUN RUN RUN RUN RUN RUN RUN RUN ] Time Iso. min 0 ] Plot cm 10 ] Range FS mW 20 ] Offset % 90 ] Pan type 1/2 1 ] Limit mW 0 Parameters ] Peak Integration ] DYN/ISO 1/2 1 ] Autolimit 0/1 1 ] Start 110 ] End 130 ] Baseline Type 8 ] Plot on 10 ] Plot Mode 101 ] DSC Integ °C 35 ] Ident. No. ] Weight mg The baseline for the integration of the curve must be specified . When there are a number of peaks on the graph, detailed evaluation of the various peaks can be carried out using the evaluation procedure. DSC INTEG °C 100 [EVALUATE] Evaluate INTEG ] Peak Integration [ [ RUN ] . Dyn/Iso 1/2 1 [ RUN ] Autolimit 0/1 1 [ RUN ] Start 110 : Start temp. for peak [ RUN ] End 130 : End temp. for peak [ RUN ] Baseline type 8 [ RUN ] Plot cm 10 [ RUN ] Plot Mode 101 [ RUN ] *** Calculating *** DSC INTEG 'C 100 -217- APPENDIX 3 Computer Program Listings NF‘O THIS FROG SUPPLIES RAD.NE.Tamb 8 8 TO SBFN3A IMPLICIT REAL‘B (A-H.O-Zi DIMENSION CP(15.15).CK(15.15).CH(15.15).A(15.15) DIMENSION ADUM(15.151,TOLD(15).DUHlIlS).DUM2(15).F(15) DIMENSION CCKIlSi.CCAPIIS).SEOLK(1$.15).R(15).D(15,1S) DIMENSION BASIlS),TEIlS),FAVEIlS).B(lS.1S) WRITEI6,1) FORMATI' ENTER RADIUS(CH), DLTH (SECS) 8 Tamb (C) ') READ(S,')RAD,DLTM,tamb IF(RAD)10,10,3 NE-l4 H-20.03 DLTM-DLTM IDTUS-l IDOEIR-lZO IDCEIM-l TAMB-tanb+273 N1-NE+1 RAD-RAD/lOO . CALL SBFNBAIRAD,N£,H,CP,CK,CH,A,B.ADUM,TOLD.DUM1,DUM2,F,CCK,CCAP, ISEOLK,R,D,BAS,TE,PAVE,N1,DLTH,IDTUS,IDDEIR.IDCEIM.TAHB) WRITE(6,4) _ [ORHAT('IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII’) WRITEI6,S) FORMAT(' DATE 15 24/5/90 ;PLOT rnoa nuru3a imanifesto”) 60 T0 9 CALL EXIT END Table A3.1 Progranntzbtunmflk '218- SUEROUTINE SBFNBAIRAD.NE.H.CP.CE.CE.A.3.ADUM,TOLO,DUM1.DUM2. 1E.CCK.CCAP.$EOLK.R.D.BAS.TE.FAVE.N1.DLTH.IDTUS,IDDEIR,IDCEIH.TAMB) C aIItIII-eeeeeetteeeseseeeeeaeeaaeaaeeeeeeeteeaesaeaeeaeeaeeaaewtt external fat.fbt.fct COMMON RI.RJ.TI.TJ.EA DIMENSION CP(N1.N1).CX(N1.N1).CH(N1.N1).A(N1.N1),B(N1,N1) DIMENSION ADUMIN1.N1),TOLDINI).DUM1(N1).DUM2(N1).P(Nl) DIMENSION CCKINI).CCAP(N1),SEOLKIN1.N1).R(N1),D(N1,N1),TE2(450) DIMENSION BAS(NE).TEINE),FAVEINE),TIMI450).!2E(4SOI.TE1(450) DIMENSION SBAIZZS).SBB(225I.SBCIZZS).FEL(1$.15),FF(1$),TPREV(1$) DOUBLE PRECISION RI.RJ,TI,TJ,Y,YY,YYY ea-96740 delcri-7.0£14 WRITEI6.')EA.DELCRI C SETTING TEE TIME/TEMP STEP CHANGE CRITERIA THEMAX-SO THEMIN-1.0 C TIME INTERVAL WILL NOT BE ALLOWED TO GO SEORTER TIMBAS VALUE TIMRAS-DLTM DUMSUM-0.0 c centhi is the max temp step allowed for centrepoint centhi-lz C NO IS COUNTER TOR XYPILES ENTRY NO-O wmmmisJJ ‘ l FORMATI' "' OUTPUT FROM SEFNBA E0. 29 "*24/S/90***') WRITEIG.‘)NE,DLTM.H C PRODUCT INFORMATION PI-3.l41$9 BRO-706.00 COND-0.0716 CAP-1547.0 c EA-78983 C......-I.IIIOIIIIIQIIIIIOIIIIIIII..-IIIIIIIIIIII'IIIIIIIf... C 4/9/85 :TINIT-ZSB FOR 'PROCESS' .303 FOR APPLE COOLING TINIT-298 C TE(NE)INITIALISED HERE FOR COOLING AS TINIT ...... CHECK UNITS TEINE)-TINIT-273 C TEINE) SET T0 0.0 FOR HEATING CURVE TEINEl-0.0 aunOLu-TINII CIIIIIIIIIIIIIIIIIIIIIIIIII!IIII.IIIIIIIIIIIIIIIIIIICIIIOIIIII NEPl-NE+1 PROGRAMME FOR GRID CALL EANGACIRAD.NE.R.NEP1) SUBROUTINE EANGACIRAD.NE.R.NEP1) EANGACIHI GRID SET UP PAYNG SPECIAL ATTENTION TO SENSITIVE AREAS DIMENSION RINEPl) WRITE (6.1) FORMATI' ENTER NE I MULTIPLE OF 7 ),RAD (CM) ') READI5.')NE.RAD SET - UP MARKERS NEPl-NE+1 ISEICE-NE/7 DELR-RAD/(B'ISEICE) IA-ISEICE IE-IA'Z IC-IA’3 IAPl-IA+1 fa 00000000000 N Table A3 . 2 Subroutine SBFNBA -219* IAP2-IA+2 ISEICZ-IA+IB+I ISEIC3-1A+IB+2 C SIMPLE GRID R(l)-0.0 DO 20 N-Z,IAP1 20 R(N)-R(N-l)+DELR DO 25 N-IAP2.ISEIC2 25 R(N)-R(N-1)+DELR/2 ‘ DO 23 N-ISEIC3.NEP1 23 R(N)-R(N-l)+DELR/4 C DELETE THE FOLLOWING TWO COMMENT FLAGS TO GET GRID PRINTOUT C DO 40 N-l,NEPl C40 WRITE(6.')N.R(N) WRITE(5.‘)NE.RAD READ(5.')NE.RAD I?(NE-7)50.2,2 DO 92 NNN-1,NEP1 WRITE(O,';NNN,R(NNN) 92 CONTINUE C DLTM IN SEC C INITIALISATION DO 30 I-1,NEP1 CCAP(I)-0.0 CCK(I)-0.0 DO 30 J-1,NEP1 CK(I,J)-0.0 CP(I,J)-0.0 CH(I.J)-0.0 SEOLK(I.J)-0.0 B(I.J)-0.0 A(I,J)-0.0 30 CONTINUE c tD0.......It.flfiitflfififiiflififittftfififi...... DO 10 N-1.NE nrxnnnn C MAKE OUT CONDUCTION MATRIX CR CCK(N)-(4'PI'COND'(R(N+l)"3-R(N)'*3))/(3'(R(N+l)-R(N))**2) C WRITE(€,')N.CCR(N) CK(N.N)-CCR(N)+CK(N.N) CK(N.N+l)-CCK(N)'(-l)+CK(N,N+l) CK(N+1,N+l)-CCK(N)+CK(N+1.N+1) c WRITE(6.42) C42 FORMAT(' NOW WRITING K MATRIX ') C DO 41 IX-leEPl C CAPACITANCE MATRIX CP CCAP(N)-PI‘REO'CAP/(lS'(R(N+l)-R(N))'*2) C WRITE(5.')N,CCAP(N) CP(N,N)-(2‘(R(N+l)"5)-20'(R(N+l)"2)*(R(N)'*3)+30*R(N+l)*(R(N)** +4)-12*(R(N)"5I)‘CCAP(N)+CP(N.N) CP(N,N+l)-(3'(R(N+1)"S)-S'(R(N+1)"4)*R(N)+5*R(N+l)*(R(N)*'4)-3* +(R(N)"5))‘CCAP(N)+CP(N.N*1) CP(N+1,N)-(3'(R(N+l)"$)-S'(R(N+1)**4)*R(N)+S*R(N+l)*(R(N)**4)-3* +(R(N)"5))‘CCAP(N)+CP(N+1,N) CP(N+1,N+1)-(12'(R(N+1)"5)-30'(R(N+1)"4)*R(N)+20*(R(N+1)**3)*(R +(N)*'2)-2'(R(N)**5))‘CCAP(N)+CP(N+1.N+1) C WRITE(6.44) C44 FORMAT(' NON SHOWING TEE CAPACITANCE MATRIX ') C DO 45 IX-1,NEP1 Table A3. 2 Subroutine SBFN3A -220- C45 NRITE(5,')(CP(IX.JX).JX-loNEPl) 10 CONTINUE c WRITE(6.46) 46 PORMAT(’ HAVE SUREACED AFTER STT. NO. 10 ') C CONVECTION MATRIX FOR T(E+1) : 3(1) CH CH(NEP1,NEPI)O4‘PI'B'(RAD"2) C FORCE TERH ANOIS AG AN DROMCRLA EP(NEP1)-4'PI‘H'TAHB'(RAD"2) C CROCH SUAS HAITRISI NA GCOTHROHOIDI CALL ARRAY‘2.NEP1.NEP1.NEP1pNEPlpSBApcx) C WRITE(51')CR(8,8),53A(54),CR(1,1),SBA(1) DO 64 IH-1,NEPI IHIR-NEP1+(IH-l)‘NEPI WRITE(6,')CR(NE?1,IH),SBA(IHIR) 4 CONTINUE CALL ARRAY(2.N1'N1.N1oNl.SBB.CH) CALL GHADD(SBA.SBB.$BC,NEP1,NEP1) CALL ARRAY(1,NEP1,NEP1,N1,NI.SEC.SEOLK) TEST-4‘PI*R‘(RAD"2)+CK(NEPI.NEPI) C WRITE(6,')TEST,SEOLR(NEP1,NEP1),CK(2,2),SEOLK(2,2) DO 99 I-1,450 153 DO 11 J-1.NEPI DO 11 JJ-IINEPI B(J,JJ)-2'CP(J.JJ)/DLTM 0‘0 C WRITE(6.')J.JJ.3(J.JJ) ll CONTINUE CALTERED eeaeeeeeeeeeeee CALL ARRAY(2,NEP1,NEP1,N1,N1,SBA,SEOLK) CALL ARRAY(2.NEP1.NEPI,N1,N1,SBB,B) CALL GMADD(SBA.SBB,SBC,NEP1.NEP1) CALL ARRAY(l,NEPl,NEPl,NI,Nl,SBC.A) c WRITE(6.60) 60 FORMAT(' BACK FROM GMADD') C CALL GMSUB(B,SEOLK,D,NEPI.NEP1) C TIME ITERATION ' WW-TE(NE) :91:-l)l!1,180,181 C NIALU 180 tan-0.0 DO 81 INC-1.NEP1 TOLD(INC)-TINIT DUMl(INC)-0.0 DUR2(INC)-0.0 DO 81 J-l.NEP1 81 CONTINUE Ctfiifififlfififittfiflfittfifltfltfitfifitflttttfittittttttifiiltitttitt'titti C H E A T G E N E R A T I O N T E R M 181 DELCRI-delCti DO 82 N-1,NE NPl-N+l RI-R(N) RJ-R(NP1) TI-TOLD(N) TJ-TOLD(NP1) c write(6.7l) CON-4'PI'DELCRI/(RJ-RI) 71 EORMAT( ’ NOW AT STT 71 ') CALL DQGBZ(RI.RJ.PCT.Y) Table A3.2 Subroutine SBFNBA '221- c WRITE(6.74) 74 FORMAT(' IST TRIP FROM DQG32 ' ) CALL DQG32(RI,RJ.FAT.YY) c WRITE(6.73) 73 FORMAT(' BACK ROM DQG32( 2ND TIME) ') CALL DQG32(RI.RJ.FBT,YYY) c WRITE(5,')Y.YY.YYY FEL(N.N)-CON‘(YYY-YY) FEL(NP1.N)-CON'(YY-Y) C NRITE(6,*)Y,YY.YYY,FEL(N,N),FEL(NP1,N)'CON 82 CONTINUE C ADDITION DO 83 N-2.NE NMl-N-l P(N)-FEL(N.NM1)+EEL(N.N) 83 CONTINUE E(l)-PEL(1.1) t(NPI)‘EEL(NP1.NE)+EE(NP1) C..............................*.........C‘.....**'.*.. 158 TPREV(J)-TOLD(J) CALL ARRAY(2.NEP1.NEP1,N1,N1,SBA.3) CALL GHPRD(SBA.TOLD,DUH1,NEP1.NEP1,1) CALL GMADD(DUM1.E.DUN2,NEP1,1) c WRITE(6.61) 61 FORMAT('BACK FROM GMPRD AND GMADD’) DO 13 JJ-1,NEP1 13 ADUH(J,JJ)-A(J.JJ) c WRITE(6.47) 47 FORMAT(' AM NOW PAST STT.13 AFTER SETING UP 3 MATRIX ') CALL ARRAY(2,NEP1.NEP1,N1,N1,SSA.ADUM) NSIMQ-O C WRITE(6,62) 62 FORMAT(' AG GABHAIL ISTEACH GO SIMQ') CALL SIMQ(SBA.DUM2,NEP1.KS) C WRITE(5.49) 49 FORMAT(' TRIP NO. TO SIMC 8 KS VALUE ARE: ') _ NSIMQ-NSIMQ+1 C WRITE(6.')NSIMQ,RS DO 14 J-1.NEP1 TOLD(J)-DUM2(J) C TOLD(J)-2'DUM2(J)-TOLD(J) DUM2(J)-TOLD(J)-273 14 DUMSUM-DUMSUM+DUM2(J) DUMSUM-DUMSUM/NEPl DUM-ABS(DUMSUM-SUMOLD) SUMOLD-DUMSUM DUMSUM-O C I? (DUM-TEEMIN)151.151,152 C151 DLTM-DLTM'Z C WRITE(6.101) C101 FORMAT(' DLTM DOUBLED,CRECEING TIMBAS ') C IF(DLTM-TIMBAS)1S4.154.159 C159 DLTM-TIMBAS C NEED TO REDO TRIS STEP ALSO c WRITE(6.102) 102 PORMAT(’ DLTM POUND TOO BIG.RESET TO TIMEAS ') C GO TO 160 cent-told(NEP1)-tprev(NEP1) Table A3.2 Subroutine SEFNBA -222- i£(cent-centhi)154.153,15 IF (DUM—TMEMAX)154,155.155 DLTM-DLTM/Z THIS TIME STEP DO 156 J-1.NEP1 TOLD(J)-TPREV(J) GO TO 153 D0 15 IL-1.NE TE(IL)-CENTERPOINT TEMPERATURE TEtIL)-(DUMZ(IL)+DUM2(IL+1))/2 TEOCHT-TAMB-273 CONTINUE TAM-DLTM/2+tan C write(6,')cent,i iftcent.gt.centhi)dltm-dltm/Z C%t%§t¥§\%t\§ DIAGNOSTIC CHECK $4884t$§8$8§84 UAIR-TAM/3600 TIME-TAM/60 - 12 CONTINUE NO-NO+1 ‘ IF (DUM2(1).GT.1000)GO TO 77 I? (DUM2(NEP1).GT.1000)GO TO 77 write(2,89)time.dum2(13) WRITE(4.89)TIME.DUM2(4) WRITE(7,89)TIME.DUM2(7) WRITE(8.89)TIME,DUM2(10) WRITE(6.')TIME.DUM2(1).DUM2(NEP1) WRITE(9,89)TIME.DUM2(Nep1) wtite(3.89)time,duml(1) TIM(NO)-TIME C SET UP MATRIX TO HOLD CENTRE TEMP...TE1(I) TEl(NO)-DUM2(1) TE2(NO)-DUM2(NEP1) ERR-WW-TE(NE) C""" HEATING SECTION ""' IF(ERR-l.0)99.95.95 C.......‘O.....fifififitfifiitiittfiifiitiiflfiflfi'..ttfiiiit"fl'.‘ nwnnnnnnn U‘ 0“.“ 040‘ ‘- mmmmm OIOMmN U 0 .e Ln 0 0000 Eeeeee COOLING SECTION eeeeeeeeee CS IF(ERR-l.0)95,99,99 93 AERR-DA3$(ERR) IF(AERR-1.0)99.95.95 C.fl......flflfltfltifitififitifi'tfififlifitii..flflfiflifiittflfififlflflfifififi 99 CONTINUE 77 WRITE(5.88) 88 PORMAT(' CONDITIONS:TINIT,TMEDIUM,ZVALUE,H,RAD,NE,TSTEP ,TAMB ') TINIT-TINIT-273 WRITE(6.')TINIT.TEOCRT.Z,R.RAD.NE.DLTM.TAMB 96 FORMAT(' UNSTABLE FOR THIS CONFIGURATION: TRY AGAIN! ') C WRITE(8,*)NO C WRITE(9,')NO C DO 84 I-loNO C WRITE(9.89)TIM(I).TE1(I) C NRITE(9.89)TIM(I).TE2(I) 80 PORMAT(IX,P7.2,',',E9.4) 89 FORMAT(1X,E12.6.'.',F12.6) C84 CONTINUE RETURN 9S WRITE(6.96) WRITE(6,*)RAD,NE,H,DLT.,TAM8,NW,DUM2(NE) fikuiUe4A3.2 Subroutine SEETEMX -223- RETURN END Table A3.2 Subroutine SBFNBA -224- FUNCTION FAT(R) COMMON RI.RJ,TI,TJ.EA DOUBLE PRECISION RI,RJ,TI,TJ,FAT TEO-((RJ-R)'TI+(R-RI)*TJ)/(RJ-RI) GAS - 8.31434 PAT- (R"3)‘(EXP(-EA/(GAS'TEO))7 RETURN END Table A3. 3 Subroutine FAT -225- ruwcrrou 937(3) cannon RI.RJ.TI.TJ,EA DOUBLE PRECISION RI,RJ,TI,TJ.FBT TEO-((RJ—R)'TI+(R-RI)*TJ)/(RJ-RI) GAS - 8.31434 PBT-(EXP(-EA/(GAS'TEO)))*RJ‘(R**2) RETURN END Table A3.4 Subroutine PST -226- FUNCTION FCT(R) COMMON RI.RJ.TI.TJ.EA DOUBLE PRECISION RI,RJ,TI,TJ.FCT TEO-((RJ-R)*TI+(R-RI)’TJ)/(RJ-RI) GAS I 8.31434 FCT-(EXP(-EA/(GAS'TEO)))*RI*(R'*2) RETURN END ThbhaAB.5 Suhnmfiinelxfl? -227- 21 100 110 120 125 130 140 SUBROUTINE ARRAY (MODE.I,J.N.M.S.D) IMPLICIT REAL'B (A-Hpo-Z) DIMENSION S(l).D(1) NI-N-I WRITE(6.21) FORMAT(' BEANNACHTAI O ARRAY ') IF(MODE-1)100.100.120 IJ-I'J+l NM-N*J+1 DO 110 K-1,J NM-NM-NI IJ-IJ-l NM-NM-l D(NM)-S(IJ) GO TO 140 IJ-O NM-O DO 130 K-1.J DO 125 L-1,I IJ-IJ+1 NM-NM+1 S(IJ)-D(NM) NM-NM+NI RETURN END Table A3.6 albrcutjne ARRAY -228- SUBROUTINE SIMO(A.B.N.KS) SUBROUTINE SIMO PURPOSE OBTAIN-SOLUTION OF A SET OF SIMULTANEOUS LINEAR EQUATIONS, AX-B USAGE CALL'SIMQ(A.B.N.K5) DESCRIPTION OF PARAMETERS A - MATRIX OF COEFFICIENTS STORED COLUMNWISE. THESE ARE DESTROYED IN TEE COMPUTATION. THE SIZE OF MATRIX A IS N BY N. B - VECTOR OF ORIGINAL CONSTANTS (LENGTH N). TEESE ARE. REPLACED BY FINAL SOLUTION VALUES, VECTOR X. N - NUMBER OF EQUATIONS AND VARIABLES. N MUST BE .GT. ONE. KS - OUTPUT DIGIT 0 FOR A NORMAL SOLUTION 1 FOR A SINGULAR SET OF EQUATIONS REMARKS MATRIX A MUST BE GENERAL. IF MATRIX IS SINGULAR , SOLUTION VALUES ARE MEANINGLESS. AN ALTERNATIVE SOLUTION MAY BE OBTAINED BY USING MATRIX INVERSION (MINV) AND MATRIX PRODUCT (GMPRD). SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED NONE METHOD METHOD OF SOLUTION IS BY ELIMINATION USING LARGEST PIVOTAL DIVISOR. EACB STAGE OF ELIMINATION CONSISTS OF INTERCRANGINC 3333 3333 NECESSARY TO AVOID DIVISION BY ZERO OR SMALL ELEMENTS. TEE FORWARD SOLUTION TO OBTAIN VARIABLE N IS DONE IN N STAGES. THE BACK SOLUTION FOR THE OTRER VARIABLES IS CALCULATED BY SU CESSIVE SUBSTITUTIONS. FINAL SOLUTION VALUES ARE DEVELOPED IN VECTOR 8. WITH VARIABLE 1 IN 3(1); VARIABLE 2 IN B(2)........., VARIABLE N IN B(N). IF NO PIVOT CAN BE FOUND EXCEEDING A TOLERANCE OF 0.0. TRE MATRIX IS CONSIDERED SINGULAR AND KS IS SET TO 1. TRIS TOLERANCE CAN BE MODIFIED BY REPLACING THE FIRST STATEMENT. nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn IMPLICIT REAL'O (A-H.O-Z) DIMENSION A(l).B(1) FORWARD SOLUTION nnn TOL-0.0 KS-O JJ--N JY-J+1 JJ-JJ+N+1 Table A3.7. Subrcmtine SD!) -229- BIGA-0.0 IT-JJ-J DO 30 I- SEARC 000 IJ-IT+I IF(DABS( 20 BIGA-A(I IMAX-I 30 CONTINUE TEST nnn IF(DABS( 35 RS-l RETURN INTER 000 40 Il-J+N*( IT-IMAX- DO 50 K- I1-II+N I2-I1+IT SAVE-A(I A(I1)-A( A(I2)-SA DIVID nun 50 A(I1)'A( SAVE-B(I B(IMAX)- B(J)-SAV ELIMI (‘00 IF(J-N) 55 IQS-N‘(J DO 65 IX IXJ-IQS+ IT-J-IX DO 60 JX IXJX-N'( JJX-IXJX 60 A(IXJX)- 65 B(IX)-B( BACK 000 70 NY-N-l IT-N‘N DO 80 J- IA- I T-C- IB-N-J ,IC-N DO 80 K! B(IB)-B( IA-IA-N Table A3.7 J,N R FOR MAXIMUM COEFFICIENT IN COLUMN BIGA)-DABS(A(IJ))) 20.30.30 J) FOR PIVOT LESS THAN TOLERANCE (SINGULAR MATRIX) CHANGE ROWS IF NECESSARY J-Z) J J,N 1) I2) VE E EQUATION BY LEADING COEFFICIENT I1)/BIGA MAX) 8(3) E/BIGA NATE NEXT VARIABLE 55.70.55 -1) -JY,N IX -JY,N JX-1)+IX +IT A(IXJX)-(A(IXJ)'A(JJX)) IX)-(B(J)‘A(IXJ)) SOLUTION 1,NY 11h, IB)-A(IA)'B(IC) Subroutine251Ft2 -230- 80 IC-IC-l RETURN END Table A3.7 Submztjne SIM) -231- BO’C.1/7/87, 32-POINT GAUSSIAN QUADRATURE SUBROUTINE DQG32 AS PER P 303. SSP MANUAL SUBROUTINE DQG32(KL.XU.FCT.Y) FCT(X) IS THE EXTERNALLY DEFINED FUNCTION XL.XU ARE THE LOWER AND UPPER LIMITS F.S. Y IS THE RESULTANT AREA UNDER THE CURVE DOUBLE PRECISION XL.XU,Y,A,B,C.FCT write(6,1) tornat(' now in dqg32' ) A-.SDO*(XU+XL) B-XU-XL C-.4986319309247407BDO*B Y-.35093050047350483D-Z*(FCT(A+C)+FCT(A-C)) C-.49280$755772634l7DO*B write(6.2) format(' now at line 14 in dqg32 ') Y-Y+.813719736545283SD-2*(FCT(A+C)+FCT(A-C)) C-.48238112779375322D0'B Y-Y+.1269603265463103OD-l’(FCT(A+C)+FCT(A-C)) C-.46745303796886984D0'B Y-Y+.17136931456510717D-1*(FCT(A+C)+FCT(A-C)) C-.44816057788302606DO*B Y-Y+.2141794901lll334OD-1'(PCT(A+C)+PCT(A-C)) C-.42468380686628499DO*B Y-Y+.25499029631188088D-1‘(FCT(A+C)+FCT(A-C)) C-.3972418979839712000'3 Y-Y+.29342046739267774D-1*(FCT(A+C)+FCT(A-C)) C-.36609105937014484DO'8 Y-Y+.329111113881809230-1*(FCT(A+C)+FCT(A-C)) C-.3315221334651076000*B Y-Y+.36172897054424253D-l'(FCT(A+C)+FCT(A-C)) C-.2938$78786203811600'B Y-Y+.390969478935351530—1'(FCT(A+C)+FCT(A-C)) C-.2$34499544661147ODO~B Y-Y+.4165596211347337BD-l*(FCT(A+C)+FCT(A-C)) C-.21067563806531767D0*B Y-Y‘ 43826046502201906r-1'fFCT(A+C!+FCT(A-Cl) C-.16593430114106382D0'B Y-Y+.45586939347881942D-1*(FCT(A+C)+PCT(A-C)) C-.11964368112605854DO*B Y-Y+.469221995404022830-1f(FCT(A+C)+FCT(A-C)) C-.722359807913982$D-1*8 Y-Y+.478193600396374300-1'(FCT(A+C)+FCT(A-C)) Co.241538328438691580-lts Y-B'(Y+.482700442573639OOD-l*(PCT(A+C)+PCT(A-C))) RETURN END nnn nnn HO Table A3.8 Subroutine 00632 -232- 10 SUBROUTINE GMADD(A.B,R,N,M) IMPLICIT REAL'B (A-H.O-Z) DIMENSION A(1).B(1).R(l) WRITE(6.54) FORMAT(' BEANNACRTAI O GMADD ') NM-N'M DO 10 I-1,NM R(I)-A(I)+B(I) RETURN END Table A3.9 Subroutine GMADD -233- 10 SUBROUTINE GMPRD(A.B.R.N,M,L) IMPLICIT REAL'B (A-H.O-Z) DIMENSION A(l).B(1).R(1) WRITE(6.64) FORMAT(' GMPRD AGB GLAOCB 1 ') IR-O IK--M DO 10 K-1,L IK-IK+M DO 10 J-1,N IR-IR+1 JI-J-N IB-IK R(IR)-0 JI-JI+N IB-IB+1 R(IR)-R(IR)+A(JI)*B(IB) RETURN END Table A3.lO Subroutine GMPRD -234- APPENDIX 4 Bibliography Abramowitz, M. and Stegun, I.A. (1972), Handbook Of Mathematical Functions. Dover Publications, Inc., New York, NY. An Bord Bainne (1988), 1987 Annual Report. An Bord Bainne/Irish Dairy Board, Dublin. Anderson, C.A. and Zienkiewicz, O.C. (1974), Spontaneous ignition : Finite element solutions for steady and transient conditions. Trans. ASME, Journal of Heat Transfer 8, 398-404. Avonmore Co-operative (1987), Personal communication, Avonmore Co-operative Ltd., Ballyragget, Co. Kilkenny, Ireland. Bartknecht, W. (1989), Dust Explosions. Springer-Verlag, New York, NY. Beever, P. (1984) , Spontaneous ignition of milk powders in a spray-drying plant. J. Soc. Dairy Technology 37, 2, 68-71. Bishop, R.C. (1981), Self-heating of materials: laboratory testing; Runaway' Reactions, Inst. Chem. Eng. Symposium Series No. 68, 2/H:1, Rugby, UK. Boddington , T., Gray, P. and Harvey, 0.1. (1972), Thermal theory of spontaneous ignition. in .bodies of arbitrary shape. Phil. Trans. Royal Society, London 270, 467-505. Boddington, T., Gray, P. and Scott, S.K. (1981), Correction of kinetic data in non-isothermal reactions with non-uniform temperatures : analytical treatments for spherical reactant masses. Proc. Royal Soc. London A 378, 27-60. Boddington, T., Gray, P. and Wake, G.C. (1977), Criteria for thermal explosions ‘with and. without reactant consumption. Proc. Royal Soc. London A 357, 403-422. Boddington, T., Gray, P. and Walker, I.K. (1980), Exothermic systems with diminishing reaction rates : -235- temperature evolution, criticality and spontaneous ignition in the sphere. Proc. Royal Soc. London A 373, 287-310. Boddington, T., Gray, P. and Walker, I.K. (1981), Runaway reactions and thermal explosion theory. Runaway Reactions, Inst. Chem. Eng. Symposium No. 68, l/C,:l, Rugby, U.K. Bowes, P.C. (1984), Self-heating: Evaluating And Controlling The Hazards. Department of the Environment, Building Research Establishment, London, England. Bowes, P.C. and Townsend, 8.8. (1962), Ignition of combustible dusts on hot surfaces . British Journal of Applied Physics 13, 105-114. Carr, R.L. (1976), Powder and granule properties and mechanics. Gas-Solids Handling in the Processing Industries, J.M. Marchello and A Gomezplata (editors), Marcel Dekker Inc., New York. ' Carrie, H.S., Walker,I.K. and Harrison, W.J.(1959), The spontaneous ignition of wool.II. A new process for wool removal from sheepskin pieces. J. Appl. Chem. 9, 608-615. Carslaw, H.S. and Jaeger, J.C. (1959), Conduction Of Heat In Solids ( 2nd Edition). Oxford University Press, Oxford, U.K. Chambre, P.L. (1952), On the solution of the Poisson- Boltzmann equation with application to the theory of thermal explosions. J. Chem. Phys. 20, 1795-1797. Clough, R. W; (1960), The finite element in plane stress analysis Proc. 2nd A.S.C.E. Conference on Electronic Computation, Pittsburgh, PA, 345-378. Crank, J. and Nicolson, P. (1947), A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Proc. Cambridge Phil. Soc. 43, 50-67. Davies, A.J. (1980), The Finite Element Method, A First Approach, Oxford University Press, Oxford, U.K. -236- Drysdale, D.D. (1985), An Introduction To Fire Dynamics. John Wiley and Sons, New York. * Drysdale, D.D. (1980), Aspects of smouldering combustion. Fire Prevention Science and Technology No. 23, Fire Protection Association, London. Duane, T.C. and Synnott, E.C. (1981), Effect of some physical properties of milk powders on Minimum Ignition Temperatures. Runaway Reactions, Inst. Chem. Eng. Symposium Series No. 68, 2/J:1, Rugby, UK. Foley, J., Buckley, D.J. and Murphy, M.F. (1974), Commercial Testing and Product Control in the Dairy Industry, Dairy Technology Department, University College, Cork, Ireland. Frank-Kamenetskii, D.A. ( 1939) , Temperature distribution in reaction vessel and stationary theory of thermal explosion. Journal of Physical Chemistry (USSR) 13, 738- 755. Frank—Kamenetskii, D.A. (1969), Diffusion .And. Heat Transfer In Chemical Kinetics ( 2nd Edition ) . Plenum Press, New York, NY. Golden Vale (1984), Personal communication, Golden Vale Co-operative Ltd., Rath Luirc, Co. Cork, Ireland. Graviner (1990), Explosion protection systems. Food Processing, January, 47. Gray, P. and Harper, M.J. (1958), The thermal theory of induction periods and ignition delays. Proceedings Seventh International Symposium of Combustion, London. Gray, P. and Lee, P.R. (1967A), Thermal Explosion Theory. In Oxidation and Combustion Reviews 2 [ ed. C.F.H. Tipper], 1-183, Elsevier, Amsterdam, The Netherlands. Gray, P. and Lee, P.R. (1967B), Proceedings Eleventh International Symposium of Combustion, Pittsburgh, PA, 1123. -237- Heldman, D.R. and Singh, R.P. (1981), Food Process Engineering ( 2nd Edition) . The AVI Publishing Company, Westport, Connecticut. Hrenikoff, A. (1941), Solution of problems in elasticity by the framework method. J. Appl. Mech. A8, 169-75. Institution of Chemical Engineers (1977), User Guide to Fire and Explosion Hazards in the Drying of Particulate Matter. I. Chem. E., Rugby, U.K. International Dairy Federation (1987), Recommendations for fire prevention in spray drying of milk powder. IDF Bulletin No. 219, IDF, Brussels, Belgium. Irish Department of Labour (1987), Health and Safety Authority, Notice of contravention of Safety Regulations on Explosion Prevention to Milk Powder’ Manufacturers, Government Press Office, Dublin, Ireland. Kordylewski, W. (1980), Critical conditions for thermal ignition of porous bodies. Combustion and Flame 38, 103- 105. Kreyszig, E. (1967), Advanced Engineering Mathematics. John Wiley And Sons, New York, NY. Larkin, J.W. (1984), Thermal diffusivity estimation from thermal process data. Ph. D. dissertation, Michigan State University, East Lansing, Michigan, U.S.A. Lovric, T., Pilizota, V. and Janekovic, A. (1987), DSC study of the thermophysical properties of aqueous liquid and semi-liquid foodstuffs at freezing temperatures, J. Food Science, 52, 3, 772-776. MacCarthy, D.A. (1983), The effective thermal conductivity of skim milk powder. Proceedings Third International Conference on Engineering in Food, Dublin, Ireland, 1, 527-538. Martin, T. (1987), Personal communication, Food Engineering Department, University College, Cork, Ireland. -238- McHenry, D. (1943), A lattice analogy for the solution of plane stress problems. J. Inst. Civ. Eng. 21, 59-82. Merzhanov, A.G., Barzykin, V.V., Abramov, V.G. and Dubovitskii, F.I. (1961), Thermal explosion in the liquid phase with heat transfer by convection only. Russ. J. Phys. Chem. 35, 1024-1027. Misra, R.N. and Young, J.H. (1979), The finite element approach for solution of transient heat transfer in a sphere. Trans. ASAE 24 , 944-949. Moshenin, N.N. (1980), Thermal Properties Of Foods And Agricultural Materials, Gordon and Breach, New York, NY. Myers, G.E. (1971), Analytical Methods In Conduction Heat Transfer. McGraw-Hill Book Company, New York, NY. O'Callaghan, D., Lafferty I., Walsh, M., Kelly, P., (1988), Explosion research at Moorepark, Co-op Ireland, June, 41-45. O'Donnell, L.P. (1988), An investigation into temperatrure and reactant consumption profiles during self-ignition of some food powders, M.Sc. dissertation, National University of Ireland, Dublin, Ireland. O'Mahony, J.G. and Synnott, E.C. (1988), Influence of sample shape and size on self-ignition of a fat-filled milk powder. J. Food Engineering 7, 271-280. Palmer, R.N. (1957), Smouldering combustion of dusts and fibrous materials. Combustion and Flame 1, 129. Quinn, J.R., Raymond, D.P. and Harwalkar, V.R. (1980), Differential Scanning Calorimetry of meat proteins as affected. by' processing treatment. J. Food. Science 45, 1146-49. Raemy, A. and Iambelet, P. (1982A), A calorimetric study of self-heating in coffee and chicory. J. Food Technology 17, 451-460. -239- Raemy, A. and Loliger, J. (1982B), Thermal behaviour of cereals studied by heat flow calorimetry. Cereal Chemistry 59, 3, 189-191. Raemy, A., Lambelet, P. and Loliger, J. (1985C), Thermal analysis and safety in relation to food processing. Thermochim. Acta 95, 441-446. Raemy, A. and Loliger, J. (1985B), Self-ignition of powders studied by high pressure Differential Analysis. Thermochim. Acta 85, 343-346. Raemy, A., Michel, F. and Iambelet, P., (1985A), Thermal behaviour of foods studied by Differential Thermal Analysis (DTA) and Differential Scanning Calorimetry (DSC). Calorimetrie et Analyse Thermique 15, 11-18. Raemy, A. and Schweizer, T.F. (1983), Thermal behaviour of carbohydrates studied by heat flow calorimetry. J. of Thermal Analysis 28, 95-108. Rice, O.K., Allen, A.O. and Campbell, H.C. (1935), The induction period in gaseous thermal explosions. J. American Chem. Soc. 57, 2212-2222. A Rothbaum, H.P. (1963), Spontaneous combustion of hay. J. Applied Chemistry 13, 291-302. Segerlind, LHE. (1976), Applied Finite Element Analysis. John Wiley and Sons, New York, NY. Semenov, N.N. (1928), Theories of combustion processes. 2. Phys. Chem. 48, 571-582. Shouman, A.R. and Donaldson, A.B. (1975), The stationary problem of thermal ignition in a reactive slab with unsymmetric boundary temperatures. Combustion and Flame 24, 203-210. Simchen, A.E. (1964), Thermal Explosions : The error in series approximations to critical temperatures. Israel J. Chem. 2, 33-34. -240- Synnott, E.C. and Duane, T.C., (1986), Fire hazards in spray-drying of milk products. In Concentration And Drying Of Foods [ ed. D.A. MacCarthy ]. Elsevier Applied Science Publishers, London. Synnott, E.C., Glynn, J.B. and Duane T.C. (1984), Minimum ignition temperatures of a self-raising flour on a standard hot-plate. J. of Food Engineering 3, 151-160. Thomas, P.H. (1958) , On the thermal conduction equation for self-heating materials with surface cooling. Transactions of the Faraday Society 54, 60-65. Thomas, P.H. (1960), Some approximations in the theory of self-heating and thermal explosion. Transactions of the Faraday Society 56, 833-839. Thomas, P.H. (1972), Self-heating and thermal ignition : A guide to its theory and application. In Ignition, Heat Release, and Noncombustibility of Materials. ASTM STP 502, American Society for Testing and Materials, 56-82. Thomas, P.H. and Bowes, P.C. (1961A), Some aspects of the self-heating and ignition of solid cellulosic materials. British Journal of Applied Physics 12, 222-229. Thomas, P.H. and Bowes, P.C. (1961B), Thermal ignition in a slab with one face at a constant high temperature. Transactions of the Faraday Society 57, 2007-2017. Tyler, B.J. and Jones, D.R. (1981), Thermal ignition in an assymetrically' heated slab : IEffects of reactant consumption. Combustion and Flame 42, 147-156. Tyler, B.J. and. ‘Wesley, T.A.B. (1965), Numerical calculations of the critical conditions in thermal explosion theory with reactant consumption. Proc. 11th Int. Symp. on Combustion, Pittsburgh 1115-1122. United States Department of Agriculture (1989), Agricultural Statistics , USDA, U . S . Government Printing Office, Washington, D.C. United States Department of Commerce (1990), Statistical -241- Abstract Of The United States, U.S. Department of Commerce, Bureau of the Census, Washington, D.C. Wake, G.C. (1980), Criticality with variable thermal conductivity, Combustion and Flame 39, 215-218. Wake, G.C. and Jackson, F.H. (1976), The heat balance in spontaneous ignition: 7. Critical parameters in special geometries for reactions of zero order. New Zealand J. Sci. 19, 23-27. Wake, G.C. and Rayner, M.E. (1973), Variational methods for nonlinear eigenvalue problems associated with thermal ignition. Journal of Differential Equations 13, 247-256. Walker, I.K. ( 1961A), The heat balance in spontaneous ignition : Part 1. The critical state. New Zealand J. Sci. 4, 309-327. Walker, I.K. (19618), The heat balance in spontaneous ignition : Part 2. The effects of superimposed Newtonian cooling. New Zealand J. Sci. 4, 328-336. Walker, I.K. (1980), The heat balance in spontaneous ignition: 10 . Linear temperature coefficient of thermal conductivity. New Zealand J. Sci. 23, 289-292. Walker, 1.1:. and Harrison, W.J. (1965A), The exothermic gaseous oxidation of scoured wool. New Zealand J. Sci. 8, 106-121. Walker, 1.x. and Harrison, W.J. (1982), The reaction between pie wool and oxygen : 2. Rate of exothermic aerial oxidation. New Zealand J. Sci. 25, 159-166. Walker, I.K., Harrison, W.J. and Hooker, C.N. (1965B), The heat balance in spontaneous ignition : Part 3. Application of ignition theory to a porous solid. New Zealand J. Sci. 8, 319-332. Walker, I.K., Harrison, W.J. and Paterson, G.E. (1968), Ignition of wool in air: Part 3-Ignition in heated room air. New Zealand J. Sc. 11, 380-393. -242- Walker, I.K., Harrison, W.J. and Read, A.J. (1967), Ignition of wool in air : Part 1-Ignition temperatures of dry wool. New Zealand J. Sci. 10, 32-51. Walker, I.K., Harrison, W.J. and Read, A.J. (1969), The heat balance in spontaneous ignition : Part 4. An equation for ignition temperature of solids. New Zealand J. Sci. 12, 302-323. Walker, 1.x. and Jackson, F.H. (1975A), The heat balance in spontaneous ignition : 5. Influence of sample shape for reactions of zero order. New Zealand J. Sci. 18, 155-172. Walker, 1.x. and Jackson, F.H. (1975B), The heat balance in spontaneous ignition : 6. Values of delta for zero- order reactions with arbitrary Biot number. New Zealand J. Sci. 18, 173-183. Walker, 1.x. and Jackson, F.H. (1977), The heat balance in spontaneous ignition : 8. Influence of ambient temperature on the critical state. New Zealand J. Sci. 20, 245-253. Walker, 1.x. and Jackson, F.H. (1978A), The heat balance in spontaneous ignition : 9. Influence of thermal conductivity for reactions of zero order. New Zealand J. Sci. 21, 519-526. Walker, 1.x. and. Jackson, F.H. (1983), Calorimetry of oxidation reactions : 5. Kinetics of exothermicity of porous solids. New Zealand J. Sci. 26, 257-275. Walker, I.K, Jackson, F.H. and Wake, G.C. (1978B), Calorimetry of oxidation reactions : 4. Significance of temperature coefficient of reaction rate. New Zealand J. SC1. 21, 537-546 Walker, I.K., Read, A.J. and Harrison, W.J. (1977), Calorimetry of oxidation reactions : 2. Measurement of a diminishing reaction rate. New Zealand J. Sci. 20, 211- 220. Walker, I.K., Wake, G.C. and Jackson, F.H. (1978), Calorimetry (of oxidation reactions : 3. An improved calorimeter equation for gaseous oxidation of solids. New Zealand J. Sci. 21, 487-495. -243- Widmann, G. (1982), Kinetic measurements on polymers. Nordic Symposium on Thermal Analysis, Helsinki, Finland. Wilson, E.L. and Nickell, R.E.(1966), Application of finite element method to heat conduction analysis. Nucl. Engng. Des. 4,1-11. Wong, H.Y. (1977), Heat Transfer For Engineers. Longman, London. Wyeth (1980), Personal communication, Wyeth (Ireland) Ltd., Askeaton, Co. Limerick, Ireland. Wyeth (1982), Personal communication, Wyeth (Ireland) Ltd., Askeaton, Co. Limerick, Ireland. Zienkiewicz, O.C. and Parekh, C.J. (1970), Transient field problems : Two-dimensional and three-dimensional analysis by isoparametric finite elements. International Journal for Numerical Methods in Engineering 2, 61-71. '244-