Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Self-heating and spontaneous combustion of wood pellets during storage Guo, Wendi 2013

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2013_spring_guo_wendi.pdf [ 4.54MB ]
Metadata
JSON: 24-1.0073583.json
JSON-LD: 24-1.0073583-ld.json
RDF/XML (Pretty): 24-1.0073583-rdf.xml
RDF/JSON: 24-1.0073583-rdf.json
Turtle: 24-1.0073583-turtle.txt
N-Triples: 24-1.0073583-rdf-ntriples.txt
Original Record: 24-1.0073583-source.json
Full Text
24-1.0073583-fulltext.txt
Citation
24-1.0073583.ris

Full Text

 Self-Heating and Spontaneous Combustion of Wood Pellets during Storage by Wendi Guo B.A.Sc., China University of Petroleum (Beijing), 2003 M.A.Sc., China University of Petroleum (Beijing), 2006  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIRMENT FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Chemical and Biological Engineering)  The University of British Columbia (Vancouver)  February 2013 © Wendi Guo, 2013 ii Abstract Self-heating of wood pellets is a major concern during long term storage. Internal temperatures rose to 57 oC in 10 days in a wood pellet silo of 21.9 m diameter in Fibreco Inc. (Vancouver, Canada) after pellets (about 20oC) were loaded into the silo. Self-heating could lead to serious accidental fires, causing enormous damage and danger to workers. In this study, the self-heating rate at different temperatures was experimentally determined, and the thermal properties were measured for wood pellets produced in British Columbia, Canada. The factors such as moisture content, pellet age and environment temperature were investigated and their impacts on the self-heating process were analyzed. Moisture content has a significant effect on effective thermal conductivity and specific heat capacity of packed pellets, but has no effect on the self-heating at the temperature range of 30 oC to 50 oC. Pellets age and environment temperature are two major factors impacting the self-heating and off-gassing process. The self-heating rate is significantly increased at higher a temperature and eventually will lead to a thermal runaway when the ambient temperature is high enough. Experimental results show that the critical ambient temperature for thermal runaway decreases as the reactor size increases. The reaction kinetics was studied at both low temperatures (30oC to 50oC) and high temperatures (100oC to 200oC) and kinetic parameters were extracted from experimental results and correlations were developed. Based on all measured properties data and kinetics data, a two-dimensional axi-symmetric self-heating model was developed to predict the self-heating process and thermal runaway in large wood pellet silo. The influences of cooling airflow rate, wall insulation, and dimension of the storage container, ambient temperature and wind condition were studied. The results show that air ventilation inside of the silo is a very effective approach for reducing self- heating and preventing thermal runaway at ambient temperatures lower than 330 K. The critical ambient temperature for a 21 m diameter silo can be as low as 36 oC in the absence of air ventilation. The current model can be used to safe guide the design and operation of large industrial wood pellets silos. iii Preface Chapter 1 is a literature review that I conducted and wrote at UBC under the supervision of Dr. C. Jim Lim, Dr. Xiaotao Bi and Dr. Shahab Sokhansanj. Chapter 3 is based on experimental work that I designed, conducted, analyzed and wrote under the supervision of Dr. C. Jim Lim, Dr. Xiaotao Bi and Dr. Shahab Sokhansanj. A version of this chapter in combination with the literature review in Chapter 1 has been published: Guo, W., Lim, CJ., Bi, X., Sokhansanj, S., Melin, ,S.. Determination of effective thermal conductivity and specific heat of wood pellets., Fuel, Available online (http://dx.doi.org/10.1016/j.fuel.2012.08.037), September 2012. Part of the experimental results has also been presented in a conference: Guo, W., Self-heating of wood pellets during storage and thermal properties measurement, 8th World Congress of Chemical Engineering, 2009, Montreal, Canada Chapter 4 is based on experimental work conducted using TAM air isothermal calorimetry in the Institute for Research in Construction (M20) of National Research Council of Canada. The experiments were conducted by Mr. Ken Trischuk. I was responsible for preparing the pellets sample, analyzing the experimental data and writing the manuscript. A version of Chapter 4 has been submitted for publication: Guo, W., Trischuk, K., Bi, X., Lim, CJ., Sokhansanj, S.. Measurements of wood pellet self-heating kinetic parameters using isothermal calorimetry. Chapter 5 is based on experimental work that I conducted at UBC and will be submitted for publication: Guo, W., Lim, CJ., Bi, X., Sokhansanj, S., Melin, ,S.. Effect of pellets age on off- gassing emission from wood pellets in storage. This work was supervised by Dr. C. Jim Lim, Dr. Xiaotao Bi and Dr. Shahab Sokhansanj. I was also responsible for experiments preparation, data collection, data analysis and manuscript preparation. Chapter 6 is based on experimental work that I conducted at UBC under the supervision of Dr. C. Jim Lim, Dr. Xiaotao Bi and Dr. Shahab Sokhansanj. I was responsible for experiments conduction, data collection, data analysis and manuscript preparation. Part of iv this chapter has been presented in a conference: Guo, W., Determination of kinetic parameters of self-ignition wood pellets during storage , AICHE Annual conference, 2011, Minneapolis, USA. Chapter 7 is based on a simulation work I conducted under the supervision of Dr. C. Jim Lim, Dr. Xiaotao Bi and Dr. Shahab Sokhansanj. Part of this chapter has been presented in a conference: Guo, W., Simulation of wood pellets self-heating in a large scale silo storage system , CSCHE Annual conference, 2012, Vancouver, Canada.  v  Table of  contents Abstract ................................................................................................................................................. ii Preface .................................................................................................................................................. iii Table of contents ................................................................................................................................. v List of tables ......................................................................................................................................... xi List of figures ..................................................................................................................................... xiii Nomenclature .................................................................................................................................. xviii Acknowledgements ........................................................................................................................... xxi 1. Introduction and literature review ................................................................................................. 1 1.1. Background................................................................................................................................ 1 1.2. Literature review ....................................................................................................................... 4 1.2.1. Chemical and physical properties of wood pellets ....................................................... 4 1.2.2. Decomposition of woody materials during storage ..................................................... 7 1.2.3. Mechanism of the decomposition .................................................................................. 9 1.2.4. Combustion ...................................................................................................................... 12 1.3. Kinetic models and parameters ............................................................................................ 17 1.3.1. Thermal properties .......................................................................................................... 17 1.3.2. Self-heating kinetic parameters ...................................................................................... 20 1.4. Research objectives and thesis layout .................................................................................. 25 1.5. Significance and unique contributions of the research ..................................................... 26 2. Experimental setup and equipment ............................................................................................ 28 vi 2.1. Measurement of thermal properties of wood pellets ........................................................ 28 2.2. Off-gassing measurement setup and equipment ................................................................ 30 2.3. Calorimetric method for measuring heat generation rate of wood pellets. ................... 32 2.4. Lab scale experimental setup for  measuring heat generation  rate of wood pellets .... 34 2.5. Material preparation and equipment .................................................................................... 37 2.6. Temperature data acquisition system ................................................................................... 39 3. Measurement of thermal properties of wood pellets ................................................................ 40 3.1. Materials ................................................................................................................................... 40 3.2. Experiments and temperature profile across the storage vessel ...................................... 42 3.3. Assumptions ............................................................................................................................ 43 3.4. Mathematic models of packed bed thermal properties ..................................................... 46 3.5. Energy balance and analysis methods ................................................................................. 47 3.6. Results and discussion ........................................................................................................... 50 3.6.1. Effect of the pellet sizes ................................................................................................. 51 3.6.2. Moisture content effect on thermal conductivity ....................................................... 52 3.6.3. Moisture content effect on specific heat capacity ...................................................... 56 3.7. Conclusions ............................................................................................................................. 57 4. Measurements of self-heating using isothermal calorimetry ................................................... 58 4.1. Materials ................................................................................................................................... 58 4.2. Self-heating kinetics and theories ......................................................................................... 60 4.2.1. Oxygen independent reactions ...................................................................................... 60 4.2.2. Oxygen dependent reactions ......................................................................................... 61 vii 4.3. Results of the heat release rate.............................................................................................. 62 4.3.1. Temperature effect on the heat release rate ................................................................ 63 4.3.2. Moisture content effect on heat release rate ............................................................... 64 4.3.3. Age effect on the heat release rate ................................................................................ 65 4.4. Heat from oxygen independent reactions ........................................................................... 67 4.5. Kinetic and heat release rate parameters of oxygen-dependent reactions ..................... 67 4.5.1. Reaction rate constant k ................................................................................................. 67 4.5.2. Heat release coefficient A1 ............................................................................................. 72 4.5.3. Reaction enthalpy ............................................................................................................ 74 4.6. Comparison of self-heating kinetic model with experimental data ................................. 75 4.7. Conclusions ............................................................................................................................. 78 5. Measurements of wood pellets off-gassing in storage .............................................................. 79 5.1. Materials ................................................................................................................................... 79 5.2. Off-gassing and emission factor ........................................................................................... 80 5.3. Results and discussion ........................................................................................................... 81 5.3.1. Effects of reactor size and head space on the gas emissions .................................... 81 5.3.2. Effect of pellets age ........................................................................................................ 91 5.3.3. Effect of test temperature .............................................................................................. 96 5.4. Conclusions ........................................................................................................................... 101 6. Determination of self-heating using lab scale experiments ................................................... 102 6.1. Materials ................................................................................................................................. 102 6.2. Thermal runaway and spontaneous combustion ............................................................. 103 viii 6.3. Crossing point and critical oven temperature .................................................................. 107 6.4. Determination of self-heating kinetic parameters ........................................................... 110 6.4.1. Frank-kamenetskii method .......................................................................................... 110 6.4.2. Crossing point method ................................................................................................. 112 6.4.3. Numerical curve fitting method .................................................................................. 114 6.5. Prediction of the critical ambient temperature ................................................................. 118 6.6. Conclusions ........................................................................................................................... 121 7. Simulation of self-heating in wood pellets silo ........................................................................ 122 7.1. Self-heating in commercial wood pellets storage silo ...................................................... 122 7.2. Computation domain and boundary condition ................................................................ 124 7.3. Governing equations and computational module ........................................................... 126 7.3.1. Governing equations ..................................................................................................... 126 7.3.2. COMSOL modules ....................................................................................................... 127 7.4. Model parameters ................................................................................................................. 127 7.4.1. Wood pellets properties ............................................................................................... 127 7.4.2. Heat generation source ................................................................................................. 128 7.4.3. Heat transfer on the outer wall surface ...................................................................... 129 7.4.4. Heat exchange between airflow and wood pellets inside of silo ............................ 133 7.4.5. Air flow velocity generated by the fans ...................................................................... 135 7.4.6. Mesh size and time step ............................................................................................... 136 7.4.7. Criterion of thermal runaway ...................................................................................... 138 7.5. Evolution of self-heating and thermal runaway ............................................................... 138 ix 7.6. Parametric study ................................................................................................................... 140 7.6.1. Thickness of the wall .................................................................................................... 140 7.6.2. Wind speed outside of the wall ................................................................................... 141 7.6.3. Dimension of the container ......................................................................................... 143 7.6.4. Heat transfer model with forced air flow .................................................................. 144 7.7. Model verification ................................................................................................................. 150 7.7.1. Comparison with industrial silos ................................................................................. 150 7.7.2. Pilot scale experiments ................................................................................................. 152 7.7.3. Lab scale experiments ................................................................................................... 155 7.8. Conclusions ........................................................................................................................... 157 8. Conclusions and recommendations .......................................................................................... 158 8.1. Conclusions ........................................................................................................................... 158 8.1.1. Measurement of thermal properties ........................................................................... 158 8.1.2. Exothermic properties of pellets ................................................................................. 158 8.1.3. Wood pellets off-gassing during storage .................................................................... 159 8.1.4. Thermal runaway and critical temperature ................................................................ 159 8.1.5. Simulation of self-heating in wood pellet silos ......................................................... 160 8.2. Recommendations for industrial practitioners ................................................................. 162 8.3. Recommendations for further work .................................................................................. 162 References ......................................................................................................................................... 164 Appendix A Model parameters ...................................................................................................... 173 Appendix B Industrial silo information ........................................................................................ 175 x Appendix C Pilot scale experiments .............................................................................................. 176 Appendix D Program ...................................................................................................................... 181 Appendix E. Self-heaitng with free convection ........................................................................... 190               xi  List of  tables Table 1.1 List of wood pellets fires accidents over last year .......................................................... 3 Table 1.2  Chemical composition of wood pellets feedstock [10] ................................................ 4 Table 1.3 Composition and properties of pine wood (chips and pellets) and prepared chars. (Adapted from Teixeira et al.[15] with permission from Elsevier) ............................................... 6 Table 1.4 Results from explosibility testing dust from white pellets and bark pellets. (Reprinted from Melin [69] with permission from the author) ................................................... 16 Table 1.5 Published data and relationships on specific heat and thermal conductivity ........... 19 Table 1.6 The critical δ values for different geometries[101, 102] .............................................. 23 Table 2.1 Compositions of standard gases for GC calibration ................................................... 30 Table 3.1 Physical properties of wood pellet samples .................................................................. 41 Table 3.2. The ratio of evaporation heat (Evap) to total energy change (ρCpΔT) ...................... 46 Table 3.3 Reproducibility of experimental data-fitted thermal conductivity and heat capacity  .............................................................................................................................................................. 50 Table 3.4 Effect of pellet size on the effective thermal properties (3 repeats) ......................... 52 Table 3.5 Thermal conductivity of pellet A obtained from experiments and predicted from model equations ................................................................................................................................. 53 Table 3.6 Thermal conductivity of pellet B obtained from experiments and predicted from literature equations ............................................................................................................................. 53 Table 3.7 Specific heat capacities of wood pellets obtained by for pellet A and pellet B ....... 56 Table 4.1 Properties of group 1 samples (as received pellets and aged pellets) ........................ 59 Table 4.2 Properties of group 2 samples (stored at room temperature 23oC) .......................... 59 Table 4.3 Properties of group 3 samples (stored under 50oC) .................................................... 59 xii Table 4.4 Reaction rate constant and activation energy of fresh samples with different moisture contents ............................................................................................................................... 68 Table 4.5 Heat release coefficient A1 data (mW/g) of Group 1 fresh pellet samples with different moisture contents. ............................................................................................................. 72 Table 5.1 Sample information and experimental conditions for off-gassing study .................. 80 Table 5.2 Summary of emission parameters for CO2 and CO with different reactor sizes and head spaces .......................................................................................................................................... 85 Table 5.3 Summary of emission parameters for CH4 and H2 with different reactor sizes ...... 85 Table 5.4 Summary of emission kinetic parameters for CO2, CO, CH4 and H2 of pellets with different ages ....................................................................................................................................... 94 Table 5.5 Summary of emission kinetic parameters for CO2, CO, CH4 and H2 of pellets at different tested temperature ............................................................................................................. 99 Table 6.1 The steady state center temperature for 3 sizes containers at different oven temperatures (Pellet E). ................................................................................................................... 108 Table 6.2 E and ∆rhA values determined from different methods. ......................................... 114 Table 6.3 E and ∆rhA values determined for different size reactors (pellet E) ...................... 115 Table 6.4 E and ∆rhA values for different types of pellets ........................................................ 116          xiii List of  figures Figure 1.1. SEM Pictures a. ground wood particle by hammer mill with 1.7 mm screen opening with x500 magnification; b. cross section of pellets with x30 magnification ( Adapted a. from Lam [16] with permission from cIRcl and b. from Lam et.al. [17] with permission from American Chemistry Society) ............................................................................... 5 Figure 1.2 Heat generation in stored three forms  of biofuels at 50 oC. (Reprinted from Rupar-Gadd [25] with permission from Växjö University Press) . ............................................... 7 Figure 1.3 Hydrolysis mechanism scheme of hemicelluloses and cellulose. (Adapted from Jacobsen et. al. [53] with permission from Springer) .................................................................... 11 Figure 1.4 Pyrolysis behavior of hemicellulose, cellulose and lignin in TGA. (Reprinted from Yang et. al. [55] with permission from Elsevier) ........................................................................... 12 Figure 1.5 Presentation of the flaming and smoldering combustion showing the respective roles of combustible volatiles and active char produced by pyrolysis under heat flux at different conditions. (Adapted from Rowell [61] with permission from American Chemical Society) ................................................................................................................................................. 13 Figure 1.6 Relationship between δ and θ0 for cylindrical vessel [102] ........................................ 23 Figure 2.1 Experimental setup for the measurement of thermal properties of wood pellets . 29 Figure 2.2 The dimension and sampling port location of the reactors (45L, 15 L, 10L and 2L)  .............................................................................................................................................................. 31 Figure 2.3 Gas chromatograph GC-14A ........................................................................................ 32 Figure 2.4 Diagram of TAM air isothermal calorimeter and its compartment of one example channel ................................................................................................................................................. 33 Figure 2.5 Experimental setup for measuring the critical oven (environment) temperature of wood pellet. ......................................................................................................................................... 35 Figure 2.6 Reactor dimensions and thermocouple locations. ...................................................... 36 Figure 2.7 Sample preparing equipments ....................................................................................... 38 Figure 2.8 Data acquisition system of the temperature measurement ....................................... 39 xiv Figure 3.1 Length distribution of pellet samples used in thermal properties experiments ..... 41 Figure 3.2 Temperature-time data recorded over a 10 hours heating period in the experimental setup of Figure 2.1...................................................................................................... 43 Figure 3.3 Typical time variation of the temperature profiles from experiments and fitted by numerical and analytical method. ..................................................................................................... 51 Figure 3.4 Comparison of experimental results (pellet A and pellet B) and predicted single pellet thermal conductivity using Eq.(3.14), parallel distribution model and serial distribution model. .................................................................................................................................................. 55 Figure 4.1 Measurement results of sample 2 (Fresh pellet sample B with 8% w.b moisture content) in closed ampoule using a TAM air isothermal calorimeter at 30, 40 and 50oC ....... 63 Figure 4.2 Measurement results of pellets B with different moisture content (sample 1, 2, 3, 9)  .............................................................................................................................................................. 65 Figure 4.3 Measurement results for group 1 samples at 50oC in closed ampoules using a TAM air isothermal calorimeter. ‘A’, ‘B’ and ‘C’ stand for pellet type A, B and C, respectively.  .............................................................................................................................................................. 66 Figure 4.4 Effect of aging on activation energy for group 2 pellets (stored at 23oC).  x1 is the days stored at 23oC before the test. ................................................................................................. 68 Figure 4.5 Comparison of experimental reaction rate constant data k(T) and Eq. (12) for group 2 pellets (stored at 23oC). x1 is the days stored at 23 oC before the test. ......................... 70 Figure 4.6 Comparison of experimental reaction rate constant data and Eq. (4.10) and Eq.(4.11) for 50oC tests of group 2 pellets (stored at 23oC) and group 3 (stored at 50oC), respectively. ......................................................................................................................................... 71 Figure 4.7 Comparison of heat release coefficient A1 data and Eq. (13) and Eq. (14) for tests at 50oC using group 2 pellets (stored at 23oC) and group 3 pellets (stored at 50oC) respectively. x1 is the days stored before the test. ......................................................................... 73 Figure 4.8 Comparison of heat release coefficient A1 data A1(T) and Eq. (15) for group 2 pellets (stored at 23oC). ...................................................................................................................... 74 xv Figure 4.9 Comparison of tested heat release rate data and predicted by Eq. Eq.(4.16) at reaction temperature 50oC for pellets conditioned at 23oC (ambient temperature). ................ 76 Figure 4.10 Comparison of tested heat release rate data and predicted by Eq.(4.16) at reaction temperature 40 oC for pellets conditioned at 23oC (ambient temperature). ............... 76 Figure 4.11 Comparison of tested heat release rate data and predicted by Eq. Eq.(4.16) at reaction temperature 30oC for pellets conditioned at 23oC (ambient temperature). ................ 77 Figure 5.1 Effects of reactor size and head space on (a) CO2 concentration and (b) emission factors for fresh pellets at 22oC. ....................................................................................................... 82 Figure 5.2 The effects of reactor size and head space on CO concentration and emission factors for fresh pellets at 22oC ........................................................................................................ 83 Figure 5.3 The effects of reactor size and head space on O2 depletion for fresh pellets at 22oC. ..................................................................................................................................................... 84 Figure 5.4 The effect of reactor size and head space on CH4 concentration and emission factors for fresh pellets at 22oC ........................................................................................................ 87 Figure 5.5 The effect of reactor size and head space on H2 concentration and emission factor for fresh pellets at 22oC ..................................................................................................................... 89 Figure 5.6 CO2 emission factor at 22 oC vs. testing days in  10L reactor of 0% and 25% head spaces for pellet samples of different ages. .................................................................................... 92 Figure 5.7 CO emission factors at 22oC vs. testing days in 10L reactor of 0% and 25% head spaces for pellet samples of different ages. .................................................................................... 92 Figure 5.8 O2 depletion at 22 oC vs. testing days in 10L reactor of 0% and 25% head spaces for pellet samples of different ages. ................................................................................................. 93 Figure 5.9 CH4 emission factor at 22 oC vs. testing days in  10L reactor of 0% and 25% head spaces for pellet samples of different ages. .................................................................................... 95 Figure 5.10 H2 emission factor at 22 oC vs. testing days in  10L reactor of 0% and 25% head spaces for pellet samples of different ages. .................................................................................... 95 Figure 5.11 Effect of  temperature  on CO2 emission factor for fresh and aged pellets ......... 97 Figure 5.12 Effect of temperature on CO emission factor for fresh and aged pellets ............ 97 xvi Figure 5.13 Effect of temperature on O2 depletion for fresh and aged pellets ........................ 98 Figure 5.14 Effect of temperature on CH4 emission factor for fresh and aged pellets. ........ 100 Figure 5.15 Effect of temperature on H2 emission factor for fresh and aged pellets ............ 100 Figure 6.1 A set of measured temperature profiles from the self-heating reactor. ................ 104 Figure 6.2 Pellets surface after self-heating experiment. ............................................................ 106 Figure 6.3 Typical temperature profiles of three different sizes of reactors with increasing the oven temperature. ............................................................................................................................ 109 Figure 6.4 Linear regression of 0ln( ) 2ln( / )C T r   vs 01/ T  for pellets A and Pellets E using FK method. ....................................................................................................................................... 112 Figure 6.5 Linear regression of ln( / ) P dT dt  vs 1/ PT   for pellets A and Pellets E using CP method. .............................................................................................................................................. 113 Figure 6.6 Temperature profiles predicted using mean activation energy value of 78.7±0.8 kJ/mol and ∆rhA value of (4.22±2.5)*10 6 kJ/(kg s) ................................................................... 117 Figure 6.7 Boundary conditions of the container (filled with wood pellets) ........................... 119 Figure 6.8 Predicted relationships between critical temperature and container radius .......... 119 Figure 7.1 Illustration of the cylindrical silo and two dimensional axi-symmetric model ..... 123 Figure 7.2 Diagram of computation domain and boundary condition for both cases with and without airflow ................................................................................................................................. 125 Figure 7.3 Diagram of the COMSOL modules used in this simulation .................................. 127 Figure 7.4 Diagram of the airflow and temperature pattern ...................................................... 129 Figure 7.5 Estimated heat transfer coeffient between the wall surface and the amibent air (temperature 288K) under free convection condition. ............................................................... 131 Figure 7.6 Estimated heat transfer coeffcient between the wall surface (305K) and the ambient air (295K) under wind condition. Note: D is the diameter of the silo ..................... 132 Figure 7.7 Estimated convective heat transfer coefficient vs. gas ventilation velocity at different temperatures ..................................................................................................................... 134 xvii Figure 7.8 Estimated airflow velocity inside of the silo (Fiberco, Inc.) ................................... 136 Figure 7.9 The residues of different stages................................................................................... 137 Figure 7.10 Four stages of self-heating. ........................................................................................ 139 Figure 7.11 The wall thickness and wall material effects on the outside and inside wall temperature. ...................................................................................................................................... 141 Figure 7.12 The temperature profile of steel wall at different wind speeds v=(0,1,2,3 ) m/s  ............................................................................................................................................................ 142 Figure 7.13 The temperature development of inner and outer walls atdifferent wind speeds  ............................................................................................................................................................ 142 Figure 7.14 Relationship between critical ambient temperature and radius of the silo ......... 143 Figure 7.15 Simulation results of airflow effect at different ambient temperatures and airflow velocities ............................................................................................................................................ 146 Figure 7.16 Solid phase temperature profiles at different airflow rates ................................... 147 Figure 7.17 Relationship between ambient temperature and airflow velocity ........................ 148 Figure 7.18 Solid phase temperature development with periodical airflow ............................ 149 Figure 7.19 Temperatures of cable C1 inside a industrial silo (80% load)............................... 151 Figure 7.20 Simulation results for airflow velocities of 0.06 m/s and 0.03 m/s. ................... 151 Figure 7.21 Temperature profiles from pilot scale reactor experimental results. ................... 153 Figure 7.22 Temperature profiles from pilot size reactor simulation results with unlimited oxygen supply ................................................................................................................................... 154 Figure 7.23 Temperature profiles from pilot size reactor simulation results with limited oxygen supply ................................................................................................................................... 154 Figure 7.24 Comparison between experimental results and simulation results in lab scale reactor ................................................................................................................................................ 156 xviii Nomenclature A1 Heat releasing rate coefficient, J/s AS Particles outer surface area, m 2/m3(packed bed) a,b,c,d Constants Bi Biot number C Gas concentration, mol/m3 CP Specific heat capacity, J/(kg K) D Diameter, m DP Particle diameter in packed bed, m E Activation energy, J/mol Evap Evaporation heat, kJ/kg fi Emission factor of off-gassing, g(gas)/kg(pellets) Gr Grashof number H Height, m hWall  Heat coefficient between the wall and air, W/(m 2 K) hloc Convective heat transfer coefficient between pellets and air flow, W/(m2 K) jH Chilton-Colburn j-factor k Reaction rate constant LC Characteristic length, m LW The thickness of the wall, m M The moisture content of pellets, in percentage wet basis, % mV Molar density, mol/m 3 Nu Nusselt number P Pressure, Pa Pr Prandtl number. Q0 Heat rate per unit length of heater, W/m xix QWall Heat flux from the silo wall to the environment, W/m 3 QG Self-heating rate from wood pellets, W/m 3 QH The convective heat from wood pellets to airflow, W/m 3 q Heat generated from pellets, J Re Reynolds number rA Reaction rate T Temperature, K t Time,s V Volume, m3 Wwt Molecule weight X Reaction conversion in fraction x1 Age of pellets (the storage time), days ∆rh Reaction enthalpy, J/kg ∆rH Reaction enthalpy, J/mol  Greek symbols α Thermal diffusivity, m2/s β Thermal Expansion coefficient, 1/K λ Thermal conductivity, W/(m K)   Porosity, volume fraction of the gas (air) ν Advection velocity of the gas phase, m/s μ Viscosity of the gas phase, Pa∙s ρ Density, kg/m3 ψ Shape factor of the particles, 0.92 for cylindrical pellet Subscript inf Ambient w Wall surface xx m Mean transfer coefficient for submerged object loc Local transfer coefficient p Bulk properties of the solid material. s Single pellets property f Fluid property (gas property) eff Effective properties. ox Oxygen dependent reaction nox Oxygen independent reaction   xxi  Acknowledgements I would like to sincerely acknowledge my supervisors, Dr. Jim Lim, Dr. Xiaotao Bi and Dr. Shahab Sokhansanj for their professional supervision, financial support, and great encouragement throughout my PhD program. I appreciate their enlightening academic instruction and patient guidance. I would also like to acknowledge my committee members, Dr. Fariborz Taghipour and Dr. Victor Lo for their constructive criticism and valuable suggestions to improve my research. I would like also to thank all people in the Biomass and Bioenergy Research Group (BBRG) for their assistance and inspiration. My thanks also go to all the faculty and staff of the Department of Chemical and Biological Engineering for their valuable advice and assistances. I would like to extend my sincere thanks to Wood Pellet Association of Canada (WPAC) for providing the pellets for my tests, Mr. Ken Trischuk at the Institute for Research in Construction (M20) of National Research Council of Canada for assisting me in using their testing facility, Mr. Staffan Melin for his invaluable suggestions, and Fiberco Inc. for their generosity of sharing the temperature data in the silo. Finally, I would express sincere gratitude to my parents for their love and encouragement throughout my life.   1  1. Introduction and literature review 1.1. Background Making pellets from wood is a well established industry in British Columbia (BC) and in North America. Wood pellets production has been increasing rapidly over the past decade, with the total production in Canada about 1.9 million metric tonnes (t) in 2011. At a bulk density of 650 to 750 kg/m³, the energy content of pellets by volume is roughly twice as much as that of round wood loaded on a truck. Pellets from BC have a typical moisture content of 4-6% and an energy density of 18 GJ/tonne compared to 8 GJ/tonne for hog fuel at 50% moisture[1]. Wood pellets are used as a direct substitute for coal and natural gas or co-fired in small and large facilities for the production of heat and electricity. However, wood pellets are hydrophilic and have to be protected from moisture during storage and handling. The pellets are sold in bulk, one-tonne jumbo bags or in small 40 lbs (18.2 kg) bags. Most of the pellets have so far been exported to Europe and are transported in large ocean vessels in volumes of 7,000 to 20,000 tonne lots. Pellets may also be transported in bulk containers, railcars or tank trucks and are the most transportable bio-fuel available today. Wood pellets are today shipped in large bulk volumes across the oceans. Before shipping wood pellets are stored in large silos typically 15 m in diameter and 20 m in height holding 2500 -3500 tonnes for a period of one to ten weeks. The accidental fires have become the biggest concern during the storage and handling, which could cause real danger to the workers, enormous economy loss, damage to the storage structure and huge air pollution. Table 1.1 lists known wood pellets fires accidents happened last year. Among those accidents, there are mainly two categories, dust explosion and spontaneous combustion. Dust explosions are caused by fast combustion of dust powders, which would happen to any powdered combustible material in an enclosed location. Dust explosions normally occur during activities such as transporting or operating, starting with an explosion and following by visible fires and smokes. Spontaneous combustion often 2 happens from inside of a piled combustible material during storage, starting with invisible combustion and following by heavy smoke, visible charring of the pellets as well as on the wall, which may or may not be followed by an explosion. Both categories are extremely dangerous and need to be studied and prevented. However, in this study, we will only focus on the self-heating and spontaneous combustion of wood pellets during the storage. Self- heating may occur in piled biomass depending upon storage conditions (temperature, moisture content and ventilation) and storage period. The heat starts to accumulate and causes temperature rise when the heat does not get dissipated quickly from the piled wood pellets. There is some evidence that spontaneous combustion may be the consequence of the self- heating process, but supporting literature is limited. This thesis is focused on studying the self-heating mechanism at both low and high temperatures and investigating conditions that may lead to the spontaneous combustion.  3 Table 1.1 List of wood pellets fires accidents over last year Location Details Category Date Armstrong Pellets plant[2] British Columbia Blast followed by visible fire spreading through the storage building. Possible dust explosion Apr.2011 Shur Fire Energy Norwich[3] wood pellets plants storage heavy smoke, no visible fire and charring wall Possible spontaneous combustion Oct, 2011 Port of Tyne [4] South Shields, United Kingdom Fire happened inside of the wood pellets silo and very deep seated.  Spontaneous combustion Nov 2011  Georgia Biomass plant [5] Germany Described as 'flash-type explosion’, which happened during operation Dust explosion Jun 2011 Tilbury power station[6, 7] United Kingdom Huge blaze broke out at power station in a storage area containing 4000 tonnes of wood pellets Dust explosion Feb 2012 Laurinburg Nature' Earth plant North Carolina[8, 9] Fire in wood pellet storage silo and second fire one week later in the same silo. Possible spontaneous combustion Mar 2012 4 1.2. Literature review 1.2.1. Chemical and physical properties of wood pellets The majority of wood pellets are produced by milling dried wood chips, bark, planer shavings or sawdust into a fine powder. The ground material is compressed and extruded through a die system, which produces a dense cylindrical shaped pellet (6 to 8 mm in diameter). Depending on the compression technology used, the woody material is exposed to temperatures around 100 °C during the compression. The stepwise thermal compression process causes the lignin to plasticize and act like glue that provides for cohesion of the pellets. Wood pellets are essentially composed of cellulose, hemicelluloses, lignin and extractives. Table 1.2 lists typical ranges of example of the chemical composition of wood pellets feedstock in British Columbia. Table 1.2  Chemical composition of wood pellets feedstock [10] Oxygenated Compounds  in % weight Cellulose 30-40 Hemicelluloses 25-30 Lignin 30-45 Extractives 3-5 Cellulose is the major chemical component of cell wall that consists of linear chains of glucose units with an average degree of polymerization (DP) at least 9,000-10,000 in wood cellulose [11]. X-ray diffraction experiments indicate that in native cellulose, the glucose units are linked together to form a crystalline structure [12]. The presence of the crystalline structure restricts the cellulose to access water and chemicals and less reactive than other components. Hemicelluloses are another group of carbohydrate components, with major compositions of (galacto) glucomanan and arabinoglucuronoxylan (refer to xylan for short) in softwood [13]; they have a degree of polymerization (DP) of 50-300 [11], much lower than cellulose. The hemicelluloses together with cellulose provide the structural components of the cell wall. Comparing to cellulose, hemicelluloses are more reactive due to the lower DP value. Lignin is a polymer consisting phenylpropane units and linked by ether or carbon- 5 carbon bonds. Lignin is an important component in wood materials. It is bound to the sugar unit in hemicelluloses by covalent linkages and generates the lignin carbohydrate complexes [14]. Lignin provides the stiffness to the cell wall and a physical barrier to enzymatic decomposition of cellulose and hemicelluloses. The molecular structures of cellulose, hemicelluloses and lignin are listed in Appendix 1A. Table 1.3 listed the chemical composition and properties of pine wood chips and pellets[15]. Wood pellets have really similar chemical composition as wood chips. In char chips and pellets, the Oxygen and Hydrogen are extremely lower than in the wood chips and pellets, while the Carbon is extremely higher. This means that the Oxygen and Hydrogen are consumed much faster than Carbon during the charring process. Manufactured from ground wood, even though wood pellets still maintain a similar chemical composition with wood, the physical structure of wood particles and wood pellets are significantly different. Figure 1.1 shows the scanning electron microscope (SEM) pictures of wood particle surface and the cross section of a wood pellet [16, 17].  Figure 1.1. SEM Pictures a. ground wood particle by hammer mill with 1.7 mm screen opening with x500 magnification; b. cross section of pellets with x30 magnification ( Adapted a. from Lam [16] with permission from cIRcl and b. from Lam et.al. [17] with permission from American Chemistry Society) In raw wood, long fibers are distributed in a parallel order (Figure 1.1a), in wood pellets (Figure 1.1b) fibers are composed of short broken fibers compressed together and 6 distributed in a random order. Due to structural differences, thermal properties that are related to physical properties are expected to be different from the solids wood. Table 1.3 Composition and properties of pine wood (chips and pellets) and prepared chars. (Adapted from Teixeira et. al.[15]. with permission from Elsevier)  Wood chips Wood pellets Char chips Char pellets Proximate analysis (wt% dry basis) Ash 0.4 0.2 1.7 1.4 Volatile matter 83.1 84.4 4 2 Fixed carbon (by difference) 16.5 15.4 94.3 96.6 Ultimate analysis (wt% dry basis) C 49.2 50.8 92.6 92.8 H 6.3 6.4 1 1.3 N <0.1 <0.1 0.2 0.3 O 45.2 43.5 3.8 3.1 S <0.2 <0.2 <0.2 >0.2 LHV (MJ kg-1 dry basis) 18.6 19.6 33.4 33.1 Char bed bulk density (g cm-3) - - 0.13 0.37 Particle skeletal density (g cm-3) - - 1.27 1.34 Particle density (g cm-3) - - 0.33 0.66 Particle porosity (-) - - 0.74 0.51 Particle average thickness (mm) - - 5.2 4.2  7 1.2.2. Decomposition of woody materials during storage 1.2.2.1. Self-heating When storing large amounts of biomass, there is a high tendency of decomposition with emissions of heat and toxic gases, depending on the storage condition [18-21]. A significant temperature increase pattern was reported recently in a 21.9 meter diameter wood pellets storage silo [22]. A maximum rate of temperature increase of 1.7 oC/h was observed and a maximum temperature of 60 oC at the center of the pile was recorded during storage. Significant temperature increases were also normally found in other biomass storage and temperature can easily rise to 60-70 oC for fresh wood [19] at which thermophilic fungi can still survive [23]. In many cases, the temperature did not go beyond 70 °C [24] and there was thus no risk of self-ignition. However, a loss of energy content in the biomass was observed.   Figure 1.2 Heat generation in stored three forms  of biofuels at 50 oC. (Reprinted from Rupar-Gadd [25] with permission from Växjö University Press) . Self-heating phenomenon was investigated on various biofuels (mainly sawdust and pellets) in the range of 20 to 60 °C using micro calorimetry [25]. The results showed that there was a peak in energy production after about 10 days for the sample stored at 50 °C during the storage period (Figure 1.2). Beyond 10 days, the heat release level decreased and then remained at a lower and constant value for up to 74 days. Constant heat production rate was also found at 50 to 90 °C for wood pellets in Sweden [26]. After several spontaneous fires, more research has been focused on the storage problems and hazard control [27, 28]. For 8 instance, blowing air through storage to cool down the bulk material can effectively eliminate self-heating [29, 30]. However available data on heat production rate for predicting self- heating, and designing air flow are quite limited in the literature and urgently needed to be investigated. 1.2.2.2. Off-gassing In addition to self-heating problem, off gassing is another serious hazard caused by decomposition of pellets during storage, which is especially dangerous for ocean transportation. Study on wood pellets off-gassing was initiated after a fatal accident occurred onboard a vessel in the Port of Helsingborg, Sweden, in November 2006, when one seaman was killed, and a stevedore was seriously injured[31]. Ocean transportation of wood pellets in confined spaces may rapidly produce lethal levels of CO and an oxygen-deficient atmosphere that may affect adjacent access space and make those spaces dangerous for people to enter [32]. UBC Biomass and Bioenergy Research Group (BBRG) have been carrying out investigations on off-gassing of BC wood pellets during storage. Five 45-liter metal containers and ten 2- liter aluminum canisters were used to study the off-gassing of different types of biomass at different temperatures [33-35]. CO, CO2 and CH4 were the major gases emitted from the stored wood pellets at the test range, with CO2 concentration being the highest and CH4 being the lowest. The concentrations of the toxic gases in the sealed space of the reactors increased over time, fast at the beginning but leveling off after a few days [32, 34, 36]. It is thus reasonable to assume that the decomposition process is a combination of an oxygen dependent reaction group and an oxygen independent reaction group. Kuang et al. [33-35] developed a first order kinetic model and predicted the off-gas emissions and concentrations in storage containers at temperatures ranging from 10 to 55 oC. As the temperature increases, both the oxygen dependent and oxygen independent reaction  rates  increase with the increase of both CO2 and CH4 concentrations [34]. Comparing to wood pellets, more VOC emissions were found in stored sawdust and wood chips, such as hexanal and pentanal compounds [21, 37, 38]. This is probably because during drying and pelletization process, 9 sawdust material is exposed to relatively high temperatures, resulting in decomposition, modification or destruction of extractives in the sawdust [39-41] 1.2.3. Mechanism of the decomposition 1.2.3.1. Metabolism and microbial activities Biomass decomposition is a complex process involving numerous biological and chemical reactions. The microbial metabolic activity is one of the important causes of self-heating and decomposition at low temperatures [42]. For instance, the fundamental exergonic reaction, oxidation of glucose, as presented below [43], releases a heat energy of 2880 kJ and 6 mole of CO2 from the oxidation of one mole of glucose. C6H12O6 +6O2→6CO2+6H2O               ΔG=-2880 kJ/mol The microorganisms produce heat and water by breaking down the sugar, and enhance the environmental factors favoring further microbial activities [19]. This leads to deterioration in biomass quality, substance losses, health hazards due to the release of toxic VOCs [33-35] and oxygen depletion[32], and potential for further self-heating of the material. Self-heating can in some cases lead to fire when heat from other sources in combination with self-heating in the biomass pile reaches the self-ignition point. Wood is a hygroscopic material which always contains some water. The water content has a significant influence on microbial activities. The study on biomass showed that the temperature in a pile depended on the initial moisture content [19]. Drying biomass will lower the microbial activity, but rewetting can revive dormant spores [44]. Lehtikangas [28] investigated the storage effect of nine assortments of wood pellets, and found  that even though the pellets are relatively dry, temperature development, moisture content increase and microbial growth were still observed at some lots during five months storage. The increased activity in stored lots caused deterioration of the pellets quality especially durability. Some researchers also showed that humidity could significantly accelerate the decomposition of biomass by releasing heat and gas, and thus the speed to self- ignition. [19, 21, 45, 46]. 10 However, the knowledge of microbial activities during pellets storage is very limited in open literature. 1.2.3.2. Chemical degradation Chemical degradation of woody materials involves a series of oxidation reactions such as hydrolysis and pyrolysis. Hydrolysis is a chemical reaction leading to the breakdown of long chain molecules, hemicelluloses, cellulose and lignin, accompanied by the uptake of water. The hydrolysis process is very complex and has been studied mainly for the purpose of producing bioethanol [47]. The hydrolysis of cellulosic biomass is generally a rate limited reaction with a very slow conversion and is normally accelerated by acid or enzymes during bioethanol production [48]. The suggested mechanisms of hemicelluloses and cellulose hydrolysis are illustrated in Figure 1.3. Two types of hemicelluloses were observed based on the reaction rate, fast-hydrolyzing and slow-hydrolyzing [49]. Having a lower degree of polymerization and consisting of branched sugar chains, hemicelluloses are much easier to be broken down than cellulose. Long chains of hemicelluloses are broken down to shorter oligomers and then further to xylose, which breaks down to other degradation products in the followed reaction [50]. Due to the crystalline structure of cellulose and the physical protection provided by hemicelluloses and lignin, the hydrolytic reactivity of cellulose is limited [51]. The amorphous portion of cellulose, non-crystalline portion, breaks down easily to generate glucoses [52, 53], which further break down to other degradation products, including the release of CO2. 11  Figure 1.3 Hydrolysis mechanism scheme of hemicelluloses and cellulose. (Adapted from Jacobsen et. al. [53] with permission from Springer) Pyrolysis becomes significant in the high temperature environment with the absence of oxygen. Oxygen attacks the biomass surface, leads to degradation of biomass, and produces heat, solid char and volatiles, mainly CO2, CO, CH4 [54, 55]. Yang et al. [55] found different pyrolysis behaviors for hemicellulose, cellulose and lignin. Figure 1.4 shows the pyrolysis behaviors of these components using the thermogravimetric analyzer (TGA). Hemicelluloses are the easiest and fastest to be pyrolyzed among the three components, with a significant reaction rate beyond 220 oC. Cellulose is the most difficult to be pyrolized due to the crystalline structure[11]. Cellulose is pyrolyzed at 315 oC. The lignin pyrolysis is different from hemicelluloses and cellulose due to the aromatic rings within the lignin, which lead to an extremely wide reaction temperature range of 100- 900 oC [55]. Fast pyrolysis is studied as an important approach of producing biofuel and charcoal [54, 56] and slow pyrolysis is the significant cause for woody biomass smoldering combustion and spontaneous ignition[45, 57, 58]. Overall, pyrolysis is a very complex thermochemical process during which biomass goes through a series of reactions [59, 60].  12  Figure 1.4 Pyrolysis behavior of hemicellulose, cellulose and lignin in TGA. (Reprinted from Yang et. al. [55] with permission from Elsevier)  1.2.4. Combustion 1.2.4.1. Smoldering and flaming combustion Instead of burn directly, cellulose materials, such as woody materials, decompose and produce a mixture of volatiles, highly reactive char and tarry components under strong heat flux[61]. Combustion of wood materials is led by a series of oxidation reactions depending on the air flow rate and the oxidation rate. Under high oxidation rate, the combustion of flammable volatiles and tarry components creates flaming combustion (Figure 1.5). At low air flow rates, the production of gas phase can be limited or diluted, which cannot maintain the gas phase combustion [62]. Under certain conditions, the oxidation of solid phase creates glowing and smoldering combustion [58]. During this stage, the solid phase temperature increases while heat and toxic gases are released from the combustion even though a flame does not occur [63]. The smoldering reaction can be invisible developing into a rapid flame spreading combustion once the air flow increases or other condition changes which might increase the oxidation rate [64, 65]. Both flame combustion and smoldering combustion are affected by the air flow rate. Smoldering propagation rate increases under higher air flow rate and the development to flaming combustion takes place much faster [62]. 13  Figure 1.5 Presentation of the flaming and smoldering combustion showing the respective roles of combustible volatiles and active char produced by pyrolysis under heat flux at different conditions. (Adapted from Rowell [61] with permission from American Chemical Society) Due to the low permeability of solid wood, the oxidation zone remains close to surface of the solid wood [63]. Wood pellets have a greater permeability and faster oxygen diffusion [30] than solids wood, which suggests that wood pellet might have a bigger chance to smolder with a faster smoldering spread rate than wood. 1.2.4.2. Spontaneous combustion Smoldering is one of the most important phenomena that may follow spontaneous combustion. Large amount of heat and flammable gases are generated while the solid phase is slowly charring when the pyrolysis reaction rate is high enough. Spontaneous ignition occurs when the smoldering develops into flame combustion. There are a number of investigations on the self-combustion of biomass and wood residues [45, 57, 66, 67]. The results showed that at temperatures above 100 °C, wood started to pyrolyze and become charred; the black residue of pyrolyzed organic material is called charcoal. If a pile of stored material is exposed to oxygen during the ongoing pyrolysis, the smoldering will turn into an overt fire. Smoldering combustion is a slow combustion, which emits heat but no light. Smoldering does not have a defined starting point. Glowing combustion produces more heat 14 and some light and flaming combustion is what we normally would think of when we think of fire. Ignition temperature is defined as the temperature at which light is emitted [68]. Investigations showed that hay reached its ignition temperature point at or near 200 °C, and the speed of combustion depended on the air temperature and relative humidity [66]. Yang [45] found that the ignition temperature of wood was around 300 °C, and the time to reach ignition increased with the increase in the density of wood. Dust from wood pellets suspended in air has been found to ignite at 450 °C in ambient conditions with limited oxygen supply, tested according to the ASTM E 1491 standard. For a stationary sample of the same dust and with ample air supply, the ignition temperature was 225 °C as tested in accordance with the Bureau of Mines RI 5624 Standard [69] Liodakis [67] showed that retardant, such as (NH4)2HPO4 (DAP) and (NH4)2SO4 (AS) reduced considerably the ignition properties of the forest species, affecting their ignition delay (combustion time) and fire point values, defined as the minimum surface temperature at which flaming combustion (ignition) occurs. 1.2.4.3. Dust and explosion Manufactured from sawdust, bulk wood pellets always contain certain amount of sawdust (fines), which also increases during handling and transportation due to the breakage of pellets [30]. The dust generated from wood pellets ranges from 500 to a few μm in diameter [69]. The dust could cause some issues during handling. For instance, dust increases the resistance to airflow causing disparate temperature distribution [30] and the fine dust, < 100μm, could be inhaled by people causing health issues [69, 70]. Among all the influences, dust explosion is the most dangerous hazard, which has been frequently happening during the handling and operating of combustible materials involving fine powders, such as coal, sulfur powders and wood dust [71-73]. The fast combustion of fine dust causes the rapid pressure increases [71] and hot gases releasing in a confined space, which lead to the dust explosion. Controlling the source of dust and oxygen is very important in preventing dust explosion. The dust limit for safe operation is provided in the MSDS of the wood pellets [70]. The ignition source is another necessary factor that 15 contributes to the dust explosion, which could be hot surfaces, mechanically generated sparks, lightning, electric currents or flames. The explosibility for wood pellets dust was tested under different condition using Canadian pellets [69, 70], with the  results listed in Table 1.4. Even though the auto-ignition temperature of dust, 450 oC, looks pretty safe, the auto-ignition temperature of dust layer is much lower, 215-225 oC. This result can also be found in other material safety data sheet (MSDS) of other dust-contained wood materials, such as, 220 oC for sawdust from untreated softwood [74] and 220 oC for sawdust from untreated hardwood [75].   16 Table 1.4 Results from explosibility testing dust from white pellets and bark pellets. (Reprinted from Melin [69] with permission from the author)  17 1.3.  Kinetic models and parameters An accurate prediction of self-heating process and self-ignition in a pellet silo is very important in order to prevent the spontaneous fires and dust explosion, as well as to minimize the quality degradation caused by pellets self-heating. Large scale experiments are usually expensive and difficult to carry out, so mathematical modeling is an appropriate approach to predict the self-heating process and self-ignition of a pellet silo. Self-heating and ignition behaviors have been studied on solid materials, such as wood and coal. Numerous theoretical studies have been carried out; mostly using combustion models originally developed for wood and coal combustion in which the prediction of the ignition condition has been incorporated [45, 57, 58, 76, 77]. Thermal properties and thermal decomposition kinetics parameters are needed in these models. Thermal properties of the materials are considered either constant [78, 79] or varying as a function of solid conversion or temperature. Kinetics parameters also have variables depending on different reaction mechanisms [80]. Simple first-order reaction kinetics is commonly used in the modeling of thermal degradation of solids fuels. Although some more advanced ignition models also take into account very complex phenomena such as the gas-phase process [81], the main objective of this study is to predict the self-heating process of the solid phase and self-ignition condition using relatively simple models. The following parameters have to be accurately investigated in order to build up the self –heating model. 1.3.1. Thermal properties Thermal properties including thermal conductivity and specific heat capacity are very important parameters in investigating self-heating and spontaneous ignition. These parameters can be measured by several well-established methods [82] and some new methods developed lately [83]. There are some thermal property data available in the literature [84]. Table 1.5 lists  a few  thermal properties data for various forms of solid wood and for agricultural materials. For solid wood, Kollman et al. [85] proposed empirical equations for estimating thermal properties as a function of density, moisture content at room temperature for wood with moisture contents from 5% to 35%. Adl-Zarrabi et al.  [83] measured thermal properties of a dry wood piece (spruce) at 20oC, 110oC and 150oC.  The results showed that the thermal conductivity is highly 18 dependent on the fiber orientation of the wood specimen, higher thermal conductivity parallel to the grain direction than perpendicular to the grain direction, and not too much on the temperature. Gupta et al.  [86] developed an empirical relation between specific heat capacity and temperature for dried softwood particles in the temperature range of 40 to 140 oC using differential scanning calorimetry (DSC) . Thermal properties of typical solid wood materials were also reported in other literatures [86-89]. The effect of moisture content on the thermal properties of alfalfa pellets was studied by Fasina et al. [46] using the line heat source method for moisture content from 7.5% to 18% w.b. (wet basis). Their method was based on a lower order approximate solution of thermal diffusion equations. Pauner et al. [90] assumed a constant thermal conductivity of 0.17 W/(m K) and a specific heat of 2.2 kJ/(kg K) for studying self-heating of biofuel pellets. Blanchard [91] at Swedish National Testing and Research Institute also mentioned a value of 0.17 W/(m K) for 6 mm wood pellets (bulk density of 603 kg/m3), however, the experimental method was not mentioned in both papers. Even though wood pellets are made of ground wood, the structure of wood particles and wood pellets are different. Due to the structure difference (Figure 1.1), the thermal conductivity of wood pellets is expected to be different from the solid wood. Therefore, instead of directly applying those reported in the literature for solid wood, it is necessary to determine the thermal properties of wood pellets and their dependence on bulk density, moisture content and temperature. 19  Table 1.5 Published data and relationships on specific heat and thermal conductivity Researchers Method and Material Thermal Conductivity,  W/(m K) and temperature range Specific Heat, kJ/(kg K) and temperature range Kollman et al. [85] General wood 0.0002 0.024dry dry   [ room temperature] 1 2(1 0.000152 1 ( ))M M     [5% <M< 35% ] , 0.0046 0.116P dryC T  [0 oC -100 oC ] Adl-Zarrabi et al. [83] TPS method, A piece of wood   0.55 W/(m K) for parallel to grain direction 0.11 W/(m K) for perpendicular to grain direction [20, 110 and 150 oC] 1.07 kJ/(kg K) for dry wood 1.38 kJ/(kg K) for M.C. 9.5% Room temperature Gupta et al. [86] Modified Fitch apparatus and DSC, Dried softwood particles L=1mm  D=6.5 mm  0.0986 W/(m K) [37 oC]  0.00546 0.524PC T  [40-140oC] Fasina et al. [46]  Alfalfa pellets M.C. 7.5 to 18.0% 0.049 0.0082eff M   [20 oC] 21.083 0.089 0.0021PC M M    [20 oC]  Peters et al. [88]  No reference or method reported 0.2 W/(m K) for Fir wood particles (4mm) 0.35 W/(m K) for Spruce wood particles CP=1.733 kJ/(kg K) CP=1.112 kJ/(kg K) Pauner et al. [90] Reference to Blomquist [84] Biofuel pellets 0.17 W/(m K) [23 oC] 2.2 kJ/(kg K) for  biofuel pellets [23 oC] *M.C.  is moisture content in percent wet mass basis;  T is  temperature, K. 20  1.3.2. Self-heating kinetic parameters As we discussed in section 1.2.2, woody materials trend to decompose and degrade due to the biological reaction and chemical reaction, causing hazard by releasing heat and gases. The decomposition is affected by various factors, including temperature, moisture content, oxygen concentration and other environment conditions. The self-heating rate is described by the kinetic parameters, activation energy (E), pre-exponential factor (A) and reaction enthalpy (∆HR), which have been experimentally measured for woody materials [26, 90, 92]. The kinetic parameters vary according to different decomposition mechanisms [80]. In the measurement of self-heating rate, two approaches can be found in the literature, micro scale experiment (calorimetric method) [93-95] and lab scale experiment. 1.3.2.1. Calorimetric method to measure self-heating rate Self-heating of bulk pellets at low temperatures is a slow process that mainly occurs in large vessels and easily influenced by the environment [6, 8]. It is difficult to study low temperature self-heating under normal laboratory conditions due to the influence of the surrounding environment. Heat conduction calorimetric methods are often used to detect small heat production and to study the thermodynamic parameters for energetic materials in the laboratory [93-97]. In the chemical fields, instruments such as differential scanning calorimeters (DSC), accelerated rates calorimeters (ARC) and thermogravimetric analyzer (TGA) are precision equipment that are normally used for measuring small heat production. The DSC usually uses small (mg) samples exposed to a rapidly changing temperature to identify a critical point where the reaction accelerates. The ARC is an adiabatic (completely insulated) calorimeter in which reactions show up as a temperature increase. Although such methods have found significant practical use, e.g. for energetic materials, they require small amount of samples [98] and don’t normally work with wood pellets due to the particle size of pellets. For measuring the tiny heat production like wood pellets self-heating, isothermal calorimetry is the proper instrument [25, 26]. Isothermal calorimetry is much more sensitive so that measurements can be conducted at lower temperatures. For example, a DSC used 21 to study reactive chemicals can have a detection limit of about 45 mW [98] for a sample size of about 1 mg. The TAM Air isothermal calorimeter has a detection limit of about 20 µW [91, 99] and a sample size of about 10 g. The specific sensitivity of the isothermal calorimeter is thus about 104 times higher than for the DSC, which means that the TAM Air can detect a reaction at certain temperature about 100 K lower than the DSC. Instead of studying heat release at, e.g. 180 °C, one can make measurements at 80 °C which is much closer to the storage temperature of wood pellets. There are only a few references in the literatures on studying heat production from stored pellets. Rupar-Gadd [25] did some tests on pellets using the isothermal calorimetry, and found that there was an increase in the heat release after 10-30 days at 50-55 °C, where after the heat release decreased to a lower constant value for the remaining days of the experiment (up to 74 days). Wadso [26] also did some test on pellets and found that at 50-90°C, close to where self-heating starts to become a problem, there was an appreciable heat production rate at the initial period. The reaction rate after 5 hours decreased because of the decrease in the oxygen pressure. The presence of metallic objects appears to have had an impact on the self-heating process where overheating and fires have started [24, 26, 100]. However, the integral kinetic parameters have not been found in the literature. 1.3.2.2. Lab scale methods to study self-heating rate The decomposition process speed increases with increasing temperature. Self-heating rate is higher at higher temperatures. Lab scale methods become more useful at higher temperatures when the self-heating rate is high. Frank-Kamenetskii (FK) method [101-103] and the crossing-point (CP) method [95] are two widely used lab scale methods to investigate the self-heating rate and kinetic parameters at high temperatures and the self-ignition condition. Both models are based on the one dimensional energy balance before the ignition  2 2P G T T C Q t x          (1.1) where Cp is the specific heat capacity of the pellets, J/kg; λ is the thermal conductivity of the pellets , W/(m K);  t is time of the self-heating reaction, s. ,QG is the self-heating rate (rate of heat generation), in J/(m3 s). 22 Semonov [102] and Frank-kamenetskii [101, 102] introduced the earliest mathematical formulations for modeling self ignition in bulk materials. Following these two models, it has been customary to assume that the reaction rate and, hence, heat generation rate varies with temperature following the Arrhenius equation:  ( ) exp( / )G rQ T h A E RT     (1.2) where, QG(T) is the self-heating rate (rate of heat generation) at temperature T, in J/(m 3 s); ∆rh is the reaction enthalpy of the self-heating reaction, J/kg; ρ is density, kg/m3; A is a pre-exponential factor; in s-1; E is the activation energy of the self-heating reaction, J/mol; R is the gas constant, 8.314 J/(mol K) ;  T is the temperature of the pellets, K. Frank-kamenetskii (1938) model was developed based on the energy balance equation for one dimension geometry. Under the steady state, the transient temperature term is eliminated (∂T/∂t=0) and energy balance equation becomes  2 2 exp( / )r T h A E RT x           (1.3) By introducing a dimensionless parameter δ, Eq. (1.3) is converted to a dimensionless form,  2 2 e         (1.4) in which the dimensionless parameter δ is defined as  0 2 2 0 E RTrhAE r e RT       (1.5) Figure 1.6 shows the relationship between the dimensionless parameter δ and the stationary solution (θ0) of Eq. (1.4). For values of δ where Eq. (1.4) has a solution, a stationary solution (θ0) is possible and the stationary temperature distribution can be obtained. For those values of δ where Eq. (1.4) has no solution, the stationary state cannot be reached and an explosion (combustion in this situation) will take place. Therefore, the maximum δ value which gives a stationary solution (δC) describes the critical condition for the self-ignition (thermal explosion). 23 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 2.5     C  Figure 1.6 Relationship between δ and θ0 for cylindrical vessel [102] The critical δ value has been found for different geometry, including cubes, spheres and cylinders [101, 102].  Experimentally the critical oven temperature T0 (the minimum oven temperature when the ignition occurs) is determined for a particular reactor size r. Rearranging Eq. (1.5), the activation energy E is obtained by plotting 0ln( ) 2ln( / )C T r   as y coordinate and 01/ T as x coordinates. F-K method can also be used as a tool to predict the critical temperature. Once the kinetic parameters (E and ΔrhA) are found , the critical temperature can be calculated using Eq. (1.5)[101, 102]. Critical δ values were determined for different geometry, including cubes, spheres and cylinders et al, part of them are listed in Table 1.6. Table 1.6 The critical δ values for different geometries[101, 102] Body δc Infinite plane slab 0.857 Infinite cylinder, radius r 2.000 Infinite square rod, side 2r 1.720 Sphere, radius r 3.333 Short cylinder, radius r, height 2r 2.844 Cube, side 2r 2.569 Regular tetrahedron, radius of insphere r, Side 2a=2√6 r 2.228 24 Crossing-point (CP) method is another widely used method to experimentally determine E[101]. Instead of a stationary solution, CP method uses a transient solution of Eq.(1.1) by assuming that when two temperatures inside of the container are identical, there will be no conductive heat transfer between these two points, and Eq. (1.1) becomes  exp( / )P r P P T C h A E RT t          (1.6) The identical temperature between this two location points is defined as the crossing point (P). After TP and / PT t   (the crossing point, and the temperature increase rate at this point) being determined experimentally, activation energy E can be calculated from the slope by plotting ln( / ) P dT dt  versus (1/ PT ). Theoretically, CP method doesn’t need as many tests as the FK method, but the no-conduction assumption will cause a systematic error. Therefore CP method by itself cannot be used as a tool to predict the critical temperature for self-ignition[90]. Scientists at the Danish Institute of Fire and Security Technology (DBI), Denmark have investigated the self-heat and spontaneous ignition of wood pellets in storage and production lines since 2005 [90, 104-106]. They used both the Frank-Kamenetskii (FK) method and the crossing-point (CP) method to analyze the experimental data generated from small-scale baskets of cubic shape. The results show that for material HP 100 and 6mm wood pellets, activation energy measured by the CP method is about 20% less than that measured by the FK method. The measured ∆HrA value found with the CP method is about 100 times smaller than when measured with FK method. This study also showed that, the FK method gives a prediction of self-heating critical side length (diameter) of 30 m at 27°C for wood pellets storage, whereas the CP method prediction gives a much smaller critical side length, 7-8 m at the same temperature. The Swedish National Testing and Research Institute also conducted extensive work on spontaneous ignition of wood pellets. A number of small and medium scale tests have been carried out and a three-dimensional finite volume CFD model was set up to simulate the self-heating in wood pellets storage using the experimental data [91]  25 The studies on determining self-heating rate of wood pellets are very limited and most of them are focused on the time to ignition condition. The propagation of the self-heating of wood pellets hasn’t been fully understood. In addition those data are based on European pellets, which are quite different from BC pellets in size and species. Therefore, it is important to develop a method to study the spontaneous ignition process for BC wood pellets, and to investigate the effects of all factors including moisture content and bulk density etc. 1.4. Research objectives and thesis layout Wood pellets are relatively new product, and limited information on self-heating and spontaneous combustion can be found in the literature. Microbial activities are usually considered as one of the biological causes of the decomposition and self-heating, even though the moisture content is quite low. The chemical reactions, including hydrolysis and pyrolysis, also contribute to the decomposition of wood pellets. When the temperature is high enough, pyrolysis dominates the degradation process and produces significant amount of heat, flammable gases, and active char. Accidents might occur due to the enormous self-heating, such as, smoldering, flaming combustion and dust explosion. The objectives of this research are to investigate self-heating rate and heat propagation before spontaneous ignition of wood pellets in storage. . Chapter 2 presents two experimental methods, calorimentric methods and lab scale method,  to measure the self-heating rate and the kinetic parameters. Chapter 3 presents the experiments designed to measure thermal properties of wood pellets, including specific heat capacity and thermal conductivity. The effects of moisture content, porosity, density and pellets on the thermal properties are discussed in this chapter. Chapter 4 describes the Isothermal calorimetry method used to measure the heat generation rate at low temperatures (30-50 oC). Kinetic parameters, including activation energy E, reaction enthalpy ∆rh, and pre-exponential factor A, obtained by regression analysis of the experimentally measured heat production rate data are given. 26 . Chapter 5 describes the experiments carried out to study the gas emission during wood pellets storage in order to gain better understanding of the degradation processes.  The effects of the age of pellets and the scale of storage container on gas emissions are presented. Chapter 6 describes a lab scale experimental setup designed to study the self-heating and spontaneous ignition at high temperature.  The results  analyzed by FK method, CP method and a numerical curve fitting method to determine kinetic parameters, such as activation energy E and Arrhenius constant ∆rhA are presented.  A self-heating model  developed based on mass and energy transfer equations for a wood pellet silo is presented in Chapter 7. This model was then used to predict the temperature profiles during the self-heating process, and to identify the thermal runaway condition and time to reach ignition. The effects of wall insulation, outside wind speed, container dimension and inside ventilation were investigated. At the end, the model was verified by comparing the simulation results with temperature profiles from industrial silo, pilot scale experiments and lab scale experiments. The last chapter, Chapter 8, summarizes all the works in this thesis. Conclusions and recommendations for further studies are  given. 1.5. Significance and unique contributions of the research Wood pellets as fuel in large power plants is a fairly new application, with limited investigation on wood pellet self-heating and spontaneous combustion behavior. Wood pellets are stored as bulk in large steel silos 4,000 – 5,000 tonnes along their supply chain. Wood pellets are transported in bulk ocean vessels that may hold up to 10,000 tonnes. The storage periods  in silos and in ships may last several months. Self-heating and spontaneous ignition are  significant safety concerns during storage and transport of wood pellets. This thesis contributes to the fundamental understanding of the self- heating process of wood pellets at low temperatures and high temperatures. The research identifies environmental conditions under which  a  thermal runaway may occur within a bulk of stored wood pellets. The improved understanding of the heat generation and self-heating process will help to control self-heating and avoid spontaneous ignition of wood pellets during transport and storage by adjusting the environmental conditions, such as ventilation air flow and optimum storage 27 dimensions. An understanding of the self-heating process can also help in understanding the mechanism of off-gassing and decomposition of wood pellets. The developed numerical model in this thesis would predict temperature development and thermal runaway conditions for a given storage container. With the predicted critical conditions for thermal runaway, the model provides specific guidelines and information for designing a storage and transportation system. The model can also be used as part of a control strategy such as ventilation scheduling and overall management of the stored wood pelletsThe following unique contributions are offered in the self-heating research field. Semi-empirical correlations are developed to predict the effective thermal conductivity and specific heat capacity, specifically for white wood pellets (Chapter 3). The heat generation rates from wood pellets at low temperatures (30oC-50oC) are quantified. The unique semi-empirical correlation is acquired from extensive laboratory experiments and the reaction kinetics are developed to predict the overall heat generation rate during wood pellets storage (Chapter 4). This study discovered that the age of pellets is a significant factor affecting both the self-heating and off-gassing process during wood pellets low temperature (30oC-50oC) storage (Chapters 4 and 5). The numerical model embraced the fundamental knowledge of the heat and mass transfer and self- heating kinetics of both low temperature and high temperature. This unique self-heating model provides a useful tool to predict temperature development during self-heating and comprehensive description of the self-heating process (Chapter 7). The primary parameters used in the numerical model were experimentally determined specifically for B.C. produced softwood wood pellets .The results from this study apply to general softwood pellets, the behavior of bark pellets and hardwood pellets might be differ. Further research on the applicability of current results to other wood types and species is recommended.  28  2. Experimental setup and equipment The thermal properties of bulk wood pellets are important parameters for studying heat transfer behavior and temperature development patterns within the stored pellets. For the purpose of investigating and predicting the self-heating phenomenon in large-scale storage silos, the thermal properties are measured in this study. The heat generation rate and kinetic parameters of self-heating due to the biological and chemical decompositions are determined using both the calorimetric method and lab scale experiment. The experimental setup and equipment are presented in this chapter. 2.1. Measurement of thermal properties of wood pellets Figure 2.1 shows the experimental setup for measuring the thermal conductivity of wood pellets. . The container is made of mild steel, 310 mm in inside diameter and 610 mm in height. An electrical rod heater of 10 mm diameter with a 39.2 ohm electrical resistance was installed at the axis of the vessel running the whole length of the cylindrical container. The electrical rod heater was connected to a constant 27.9 volts DC power source. Four K-type thermocouples (T0, T1, T2 and T3) of 3 mm in diameter and 300 mm in length were placed inside the container with the tips of thermocouples located at around the middle of the container height.  T1, T2 and T3 were inserted from the top at 66 mm, 69 mm and 113 mm away from the axis, respectively. T0 was inserted from the side wall for monitoring the surface temperature of the heater. These thermocouples were connected to a data acquisition system to record the temperature every minute. The outside wall of the container was insulated by fiber glass insulation material to reduce heat losses. The heater and thermocouples were first inserted and fixed in the test container. A randomly packed bed was archived in this test by pouring pellets into the container slowly from the top opening until the container is fully filled. The top of the pellet pile was then leveled and the container’s cover was bolted to the top. The weight of pellets was recorded to calculate the filled volume of the container and the bulk density of loaded pellets. A data acquisition system (detail in section 2.6) was used to record temperature data. All pertinent data including bulk density, temperatures, and moisture 29 contents were stored in a computer file. When the wire and heater were all connected, heater was turned on while temperatures were recorded at an interval of 60s. The experimental results of a period of 150 minutes (9000 s) were used for data fitting to find thermal properties. The pellets were removed by tipping the container after turning the heater switch to its off position. 610 mm Data logging Electric heater controller Computer T1 310 mm T2 T3 E, Electric  heater 66mm 69 mm 113mm T0 Fibre glass insulation Fibre glass insulation  Figure 2.1 Experimental setup for the measurement of thermal properties of wood pellets  30  2.2. Off-gassing measurement setup and equipment   The gas emission of the wood pellets in storage was investigated using 4 sizes of storage containers (45 L, 15 L, 10 L and 2 L). The dimensions of the containers are listed in Figure 2.2. The 45 L container was made of mild steel, 310 mm in inside diameter and 610 mm in height  with two sampling ports. The upper port is 160 mm from the top and the lower port is 135 mm from the bottom. The 15 L and 10 L containers are made of plastic with one sampling port on the top. The 2 L container is made of glass with one sampling port on the top. All containers were sealed to air tight during the experiments. 25 ml of gas was sampled every day from each container, except for the 2 L container for which 15 ml was collected. The gas composition was analyzed by a gas chromatograph, GC-14A (Shimadzu Corporation, Japan), as shown in Figure 2.3. The GC-14A gas chromatograph is configured with 1 sample loading loop, an automated valve system for routing samples, a methanizer, one Flame Ionization Detector (FID) and one Thermal Conductivity Detector (TCD). The TCD detector provides a qualitative and quantitative analysis of H2, O2, N2 while the FID provides qualitative and quantitative analysis of CO, CO2, CH4 and lighter hydrocarbons (>C2). The GC system was calibrated weekly using three calibration gases (Table 2.1) to assure the accuracy of the concentration. Table 2.1 Compositions of standard gases for GC calibration Standard gas 1 Standard gas 2 Standard gas 3 Components Conc. Components Conc. Components Conc. CO2 6% CO2 0.50% CO2 0.10% CO 2.50% CO 0.10% CO 0.05% CH4 1.50% H2 1.00% H2 0.50% He 2% CH4 0.10% CH4 0.50% N2 3.00% Air balance Air balance O2 10%  Argon balance  31 610 mm 310 mm 45 L  Steel cylinder Rubber Gasket Sampling Port 1 Sampling Port 2 1 3 5  m m 1 6 0  m m  285 mm 255 mm 2 7 5  m m 15 L Plastic bucket Sampling Port  a. 45 L  b. 15 L  240 mm 203 mm 2 7 0  m m 10 L Plastic bucket Sampling Port  105 mm 2 6 5  m m Sampling Port 2 L glass cylinder  c. 10 L d. 2 L Figure 2.2 The dimension and sampling port location of the reactors (45L, 15 L, 10L and 2L) 32  Figure 2.3 Gas chromatograph GC-14A 2.3. Calorimetric method for measuring heat generation rate of wood pellets. Thermal Activity Monitor (TAM) of the air isothermal calorimeter (Model 3114, TA Instruments, 109 Lukens Drive, New Castle, DE 19720) was used to measure the heat generation rate of wood pellets at different constant temperatures. TAM air is an eight channel micro calorimeter from TA Instruments designed for sensitive and stable heat flow measurements and is ideal for large-scale calorimetric experiments requiring sensitivity in the milliwatt range, like self-heating of wood pellets. In the current study, measurements were conducted at 30, 40 and 50 oC, respectively. The thermostat employs circulating air and an advanced regulating system to keep the compartment temperature very stable (within ± 0.02 oC). Figure 2.4 shows the diagram of the TAM and its compartment channels). All calorimetric channels are of twin type, consisting of a sample (wood pellets) and a reference (glass bean) vessel, each with a 20 ml volume. The sample vessel can hold about 5 pellets, while the reference vessel contains glass beads of the same heat capacity as the sample. The thermostat compartment is kept at the stable temperature while the heat flux to both vessels is recorded. The heat generation rate of the wood pellets is calculated as the heat flux difference of the two vessels. The high accuracy and stability of the thermostat makes the 33 calorimeter well-suited for heat flow measurements over extended periods of time, e.g. weeks. A data set of specific heat flow (µW/g) vs time can be obtained from this equipment.  Figure 2.4 Diagram of TAM air isothermal calorimeter and its compartment of one example channel In preparation for a test, the TAM instrument was calibrated for the specific temperature and output ranges following the TAM air isothermal calorimeter calibration procedure [107, 108]. About 5 g of pellet samples were weighted and loaded into closed ampoules (with a volume of 23.5 ml) at ambient pressure, and the heat capacity of each loaded sample was calculated for each channel based on loaded mass and the Cp values. The required amount of glass beads for each reference channel was then calculated to match the heat capacity of each sample material. (The glass beads heat capacity is 0.90 kJ/(kg.K)). The required amount of glass beads was then loaded into the reference ampoule for each reference channel. After the system reached equilibrium, the sample ampoules were then placed into the calorimeter and the heat release rate of sample pellets was recorded by the Picolog recorder program in the computer by comparing heat flux of the sample ampoules and the corresponding reference ampoules. 34 2.4. Lab scale experimental setup for  measuring heat generation  rate of wood pellets Lab scale experiments were carried out using a laboratory oven (Thermal scientific, 1.7 cu. ft) to investigate the self-heating at high temperature and the spontaneous combustion phenomenon. Figure 2.5 shows the experimental setup. An open-top short cylindrical  steel container as the reactor was placed inside an oven and two thermocouples were placed inside the container , one in the middle (T1) and another one half way from the center (T2) to the wall (radial). The ends of thermocouples were kept at the same axial level from the bottom. The containers were filled with wood pellets (weight was recorded) and the oven was maintained at a constant temperature between 120 and 200 oC. The temperatures were recorded continuously by a data acquisition system (detail in section 2.6). Three sizes of steel cylindrical containers were used in this experiment, small (H 114 mm, D 105 mm); medium (H 188 mm, D 164 mm); and large (H 224 mm, D 216 mm).  Experiments were carried out to obtain the critical oven temperatures and the oven preset temperature where the temperature inside of the bulk pellets start to run away, for each size of the container by a number of trial-and-error iterations according to the FK method. First, the container filled with fresh wood pellets was put in the oven at a set temperature for at least 24 hours, and the wood pellet temperatures were recorded. If the wood pellet temperature increased continuously over the oven temperature (temperature runs away), the experiment was repeated at a lower oven temperature, otherwise, the experiment was repeated at a higher oven temperature. The minimum oven temperature at which the wood pellet temperature increases continuously to higher than the oven temperature is defined as the critical environment (or oven) temperature for the particular size of the reactor. Critical temperatures of three reactors were determined to estimate a value of the activation energy (E) by linear regression.  An self-heating reactor, a top open cylindrical aluminum reactor of 417 mm in height and 356 mm in diameter, was used for investigating the self-heating and spontaneous combustion phenomena and its critical ambient temperature. The temperature profiles were monitored using four thermocouples. The axial distances between the thermocouples and the center axis were 116 mm (T1), 58 mm (T2), 0 mm (T3) and 88 mm (T4), respectively, and their tips were maintained at the same axial level of 200 mm above the bottom. 25 kg pellets were loaded and the critical oven 35 temperature was measured following the trial-and-error method. The reactor dimensions and thermocouple locations are illustrated in Figure 2.6.  oven Wood pellets container Thermocouple 1 Thermocouple 2 Data Logging interface  Figure 2.5 Experimental setup for measuring the critical oven (environment) temperature of wood pellet.  36 2 2 4  m m 216 mm 1 8 8  m m 164mm 1 1 4  m m 105 mm 4 1 7  m m 356 mm 88 mm 58mm 116mm T1T2T3T4 Self-heating reactor Large reactor Medium reactor Small reactor 5.6 mm T1 T2 T1 T1 T2 T2 40 mm 28 mm 1 0 1  m m 9 4  m m 5 7  m m 2 0 0  m m  Figure 2.6 Reactor dimensions and thermocouple locations. 37 2.5. Material preparation and equipment Properties of wood pellets produced from different wood pellet plants might be slightly different. In order to maintain the consistency and accuracy of the experimental results, the physical properties such as moisture content, bulk density and size distribution are measured and recorded before each test. Bulk density was measured following the American Society of Agriculture and Biological Engineering (ASABE) standard S269.4 [109]. Moisture content of wood pellet was measured following the ASABE standard S358 [110] using a THELCO laboratory oven (Figure 2.7a). Wood pellets often contain small amount of moisture as manufactured, typically 4-8% wet based (w.b.). Due to the hygroscopic characteristics of wood material, even small amount of moisture content has an important influence on their physical and thermal properties. In order to investigate the moisture effect, wood pellets were conditioned to prepare pellets at different moisture contents in the lab. Pellets with lower moisture content were obtained by oven drying (THELCO laboratory oven) and pellets with higher moisture contents were obtained by moisturizing using an ESPEC LHU-113 temperature and humidity cabinet (Figure 2.7b). Besides the oven drying techniques, a freeze dryer was also used to dry the pellets while investigating the temperature effect on self-heating rate at low temperatures. The pellet size distribution was determined by measuring the diameter and length of 200 pellets from each lot using a micrometer. A Gilson sieving device model Test-Master shaking sieving device (Figure 2.7c) was used to fractionate the sample to two fractions of larger and smaller than 6.7 mm. 38   a. THELCO Precision laboratory oven (Thermo Fisher Scientific Inc., 308 Ridgefield Court Asheville, NC 28806) b. LHU-113 temperature and humidity cabinet (ESPEC North America, Inc., 4141 Central Parkway Hudsonville, MI 49426)  c. Test-Master shaking sieving device (Gilson Company Inc., P.O. Box 200 Lewis Center, Ohio 43035-200) Figure 2.7 Sample preparing equipments 39 2.6. Temperature data acquisition system The temperature data acquisition system (Figure 2.8) consists of thermocouples, USB-TC 8 channel thermocouple input module (Measurement computing Corporation, 10 Commerce Way, Norton, MA 02766, USA) and a computer. The thermocouples used in this research are K type with grounded tip and ceramic insulated lead wire, which connects to the USB-TC device. A Visual Basic program was developed to convert the input signal into temperature data and record them at preset periodic intervals.  Figure 2.8 Data acquisition system of the temperature measurement 40  3. Measurement of thermal properties of wood pellets1 Thermal properties of wood pellets, such as specific heat capacity and thermic conductivity, are important factors that affecting self-heating and spontaneous ignition of wood pellets in storage. This chapter describes experiments designed and conducted to determine the thermal properties of wood pellets. The effects of    moisture content, porosity, density of wood pellets, and pellets from different manufacturers are presented. . 3.1. Materials The effective thermal conductivity and specific heat capacity of two softwood pellet samples produced by two pellet plants in British Columbia were measured in this research.  Table 3.1 lists the physical properties of these two pellet samples. Both pellet A and pellet B were made of whitewood mostly pine. Pellet A was manufactured by Princeton Co-Generation Crop. (Princeton, British Columbia, Canada, V0X 1W0). Pellet B was manufactured by Premium Pellet Ltd. (Vanderhoof, British Columbia Canada, V0J 3A0). Pellet A was from a pellet lot with as a received moisture content of 4.6 %. Pellet B was from two lots; B1 with a as received moisture content of 4.6 % and B2 of 8.8%. The moisture content of pellets as received was measured following the ASABE standard S358 [111] using a THELCO laboratory oven. Pellets samples (A and B) with moisture content other than as received were obtained by conditioning in a humidity chamber ESPEC LHU- 113 (high moisture content) or oven drying by the THELCO laboratory oven (low moisture). The pellet size distribution of each lot was determined by measuring the lengths and diameters of 200 pellets using a micrometer. Figure 3.1 shows  the size distribution of these pellets.    1 A version of this chapter were published. Guo, W., Lim, CJ., Bi, X., Sokhansanj, S., Melin, ,S.. Determination of effective thermal conductivity and specific heat of wood pellets., Fuel, 2013. 103: P. 347-355 41  2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 0 5 10 15 20 Pellet Length (mm) F re q u e n c y  ( % )  Size Distribution for Sample 1 (Normal)  a. Pellets A 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 Pellet Length (mm) F re q u e n c y  ( % )  Size Distribution for Sample 2A (Moisture content 4.6%)  b. Original pellets B1  (Moisture content 4.6%) 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 Pellet Length (mm) F re q u e n c y  ( % )  Size Distribution for Sample 2B (Moisture content 8.8%)  c. Original pellets B2 (Moisture Content 8.8%) Figure 3.1 Length distribution of pellet samples used in thermal properties experiments  Table 3.1 Physical properties of wood pellet samples 42 Samples Pellet A Pellet B Pellet B1 Pellet B2 Density, kg/m3 640-670 690 590-620 Specific density, kg/m3 1159 1186 1168 Average diameter, mm 6 6 6 Average length, mm 10 13.1 13.6 Moisture content (as received), wt %, wet based. 4.6% 4.6% 8.8% Lc, mm 1.17 1.22 1.23 Biot Number 0.023-0.090 0.024 -0.094 0.024-0.094 Note: Lc is the characteristic length of a single pellet, calculated using Eq. (3.2); Biot number is calculated using Eq. (3.1) . 3.2. Experiments and temperature profile across the storage vessel The experimental setup for the thermal properties measurement was described in section 2.1. The test vessel is a steel cylindrical container of 310 mm in diameter and 610 mm in height with an electric heater of 10 mm in diameter and 610 mm in height installed at the axis of the vessel. After the pellets are loaded in the vessel, the electric heater is turned on and slowly heats up the pellets from the center with a heating rate of 32.56 W/m. During the heating, temperature variations with time were obtained using four K-type thermocouples, T0, T1, T2 and T3 located at 0, 66 mm, 69 mm and 113 mm away from the axis, respectively. Figure 3.2 shows a typical experimental result on the time variations of the temperature during the vessel heating over 10 hours using wood pellets A. Temperature profiles of T0, T1, T2 and T3, measured at the center, 66 mm, 69 mm and 113 mm away from the axis, respectively. Initially the entire bulk of wood pellets is at a uniform room temperature of 25 oC. The temperature at the container center (T0) where the heater was located increased immediately following the flow of current to the heater and reached the highest temperature of 75 oC after 10 hours. Temperatures at points T1, T2 and T3 remained unchanged for a period of almost 2500 s. At 9000 s, the center temperature (T0) rose to 55oC while T1, T2 and T3 lagged behind at 33oC, 31oC and 27oC respectively. The temperature profiles are used to determine the thermal conductivity and specific heat capacity by applying the energy balance across the vessel. However, the temperature profiles can be 43 influenced by many other factors, such as the heat loss from the container wall and the moisture evaporation inside of the bulk pellets. In order to eliminate those influences, only the temperature profiles of the first 9000 s were used to determine the thermal properties where the moisture evaporation heat loss and the heat loss from the container wall are neglected. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 4 20 30 40 50 60 70 80 Time(s) T e m p e ra tu re (o C )  Figure 3.2 Temperature-time data recorded over a 10 hours heating period in the experimental setup of Figure 2.1.  3.3. Assumptions The energy balance across the wood pellet container was used to extract the thermal properties from the temperature profile data. The following assumptions were  made to simplify the energy balance. a. The bulk wood pellets are treated as a continuous medium and the temperature distribution inside of pellets is uniform. The Biot number is a dimensionless number indicating the ratio of heat transfer resistance between the inside and surface of a solid in a heterogeneous heat transfer system. It is only reasonable to T0 T1 T2 T3 44 assume the system is a continuous heat transfer medium when the Biot number is smaller than 0.1 [112]. For this study the Biot number of individual wood pellets is defined as follows,  C s hL Bi    (3.1) where h is the convection heat transfer coefficient between single pellet surface and its surrounding air. h varies from 5-20 W/m2 K for wood and air when the air velocity ranges from 1.0 to 5.0 m/s [113].  λs is the thermal conductivity of a single pellet, which is typically 0.22 W/(m K) based on experimental results. LC is the characteristic length of a single pellet, equal to the volume (V) divided by surface area (AS) defined as follows,  2 2 ( / 2) 2 ( / 2) 4 2 C S V D L DL L A DL D L D          (3.2) where D is the pellet diameter and L is the pellet length. When Bi≤0.1, the variation of temperature within a single pellet is considered small, and therefore temperature distribution inside the bulk of pellets can be reasonably assumed to be uniform, as in a continuous system [112]. Biot numbers of two samples were calculated and shown in Table 3.1 to be less than 0.1. Therefore, for airflow smaller than 5.0 m/s, pellets storage system can be reasonably approximated as a continuum without temperature gradients within single pellets. Since there is no airflow within the bulk of pellets in this study, heat transfer coefficient h would be much smaller than the assumed 5 W/m2 yielding even a smaller Biot number. As a continuum, the entire bulk of pellets behaves as a thermal conductive medium. b. The evaporation of water moisture of pellets can be neglected for the initial first 9,000 s heating period, The energy balance across the bulk pellets in the vessel including the heat of evaporation of water moisture can be expressed as  2 2 1 ( ) vapPb eff T T T e t r r r C            (3.3) 45 where the left item of the equation is the energy changing rate kJ/(m3 s) while evap is the evaporation energy kJ/(m3 s). Figure 3.2 shows that within the fitting period (0 to 9000 s), the center temperature (T0) reached 55oC while T1, T2 and T3 reached 33, 31 and 26oC respectively. After the experiments, the inside surface of the container and top cap was checked for liquid water during the unload pellets. When the experiments end at 9000s, no liquid water was observed on the inside container wall nor the inside surface of the top cap but few liquid water droplet was observed on the inside surface of the top cap when the experiment end at 11 hours. For 9000s heating, no liquid water being observed inside of the sealed container after experiment so that it is reasonable to assume that the moisture in the air was equal to or less than saturated. Since no condensation was observed, we can assume that the moisture loss from the pellets is equal to the moisture gain in the air and the air inside the container is always saturated with water vapor during the heating process. Assuming the air is well mixed in a closed container, the molar density of moisture (mv , mol/m 3) in the air can be calculated from the water saturated vapor pressure using the idea gas law.  ( )sat v P RHn m V RT    (3.4) where RH is the relative humidity and Psat is the saturated vapor pressure; V is the air volume in the storage container. The evaporation heat loss Evap (kJ/m 3) over a time range of Δt can be calculated by multiplying the moisture loss and the vaporization enthalpy hfg (kJ/kg).  vap vap v wt fgE e t m W h       (3.5) where Wwt is the molecule weight of water; mv is the molar  density of the evaporated moisture. Assuming that at time zero, the vessel temperature is 25oC and relative humidity (RH) of air is 40% and, the air is saturated (RH=1) after 9000s, the evaporation energy loss Evap at various air temperatures were calculated and listed in Table 3.2. When the air temperature is 40oC (average temperature after 9000 s), the ratio of evaporation heat loss (Evap) to total energy change (ρCpΔT) is 0.0073 and the maximum ratio is 0.0082 when the temperature is at the maximum of 55oC in the fitting range (9000 s).  By neglecting the evaporation heat removal (Evap) in Eq.(3.3), a maximum error of 0.8% is introduced into the system. Therefore, it is safe to ignore the evaporation heat removal (Evap) in this study. 46 Table 3.2. The ratio of evaporation heat (Evap) to total energy change (ρCpΔT) Tair, oC Psat*, Pa RH mv, mol/m3 hfg* , kJ/kg Evap, kJ/m3 ρCpΔT Eevp/ρCpΔT 25 3169 0.4 0.511 2442 40 7384 1 2.837 2407 100.77 13650 0.00733 50 12350 1 4.598 2383 175.32 22750 0.00769 55 15760 1 5.779 2371 224.81 27300 0.00824 *Note: mv is the molar density of moisture;  the saturation pressure of water Psat and the evaporation enthalpy hfg are read directly from literature  [112].  3.4. Mathematic models of packed bed thermal properties As a well-accepted approximation, the packed bed of pellets can be simplified by assuming it is a continuous homogeneous porous system with effective thermal properties of λeff and CP. For a porous material that consists of solids and gas, the effective heat capacity is approximated by,  , ,(1 )b P s P s f P fC C C       (3.6) where CP,s,CP,f and ρs, ρf are specific heat  and density of the solid and the gas, respectively;  is the volume fraction of gas (porosity) in the system. Cp is specific heat of the bulk (gas plus solid), ρb is the density of the bulk. The density of air (ρf=1.127 kg/m 3) at 20 oC is much smaller than the density of solid (ρs=1159 kg/m 3), so that the effective heat capacity of bulk pellets does not differ from the specific heat capacity of a single pellet. The effective thermal conductivity of solid- fluid packed bed system is more complex. Thermal conductivity depends upon the porosity () and orientation of the solid particles distribution. Several theoretical models and empirical correlations are available in the literature to predict the effective thermal conductivity of the particle-gas packed bed systems [114-118]. The maximum thermal conductivity of solid-fluid system is obtained for parallel distribution when the solid layers parallel to the direction of the heat flow; and the minimum value is obtained for serial distribution when the solid layers are perpendicular to the direction of the heat flow [115, 117, 118]. 47 for parallel distribution                         (1 )eff f s       (3.7) and for serial distribution 1 1 eff f s          (3.8) where λs and λf are the thermal conductivity of solid and gas and  is the porosity; λeff is the effective bulk thermal conductivity. 3.5. Energy balance and analysis methods With the above assumptions a and b, a packed bed of pellets can be considered as a homogeneous materials with effective thermal properties of λeff and CP. For a bulk material that contains moisture, the effective heat capacity can be approximated by,  , ,(1- /100) ( /100)P P dry P MC M C M C   (3.9) where CP,dry,, and CP,M are specific heat capacity of dry material and the moisture, respectively; M is the moisture content in percentage, wet based. Cp is specific heat capacity of the bulk material (wet based). Unlike specific heat, the moisture effect on thermal conductivity is more complicated and cannot be simply represented as a mixture formula similar to Eq. (3.9). The moisture might migrate inside of the material due to the temperature difference and speed up the heat transfer process. Hence, the effective thermal conductivity appears larger with the existence of moisture, even for low moisture contents materials, such as wood pellets (4-8% for commercial pellets). With the assumption that the end effects can be neglected for a cylindrical geometry and the material is a continuous medium, the temperature distribution inside the cylindrical container during heating can be represented by the following partial differential equation with the stated boundary and initial conditions.  2 2 1 ( )b P eff T T T t r r r C           (3.10) 48 where T is temperature, K; r is radial coordinate, m; t is time, s; λeff is the effective thermal conductivity, W/(m K); ρb is the bulk density, kg/m 3; CP is the bulk specific heat capacity, kJ/(kg K); R is the radius of the vessel, m. Eq. (3.10) assumes the bulk pellet is a continuum of biomass solids and interstitial air mixture. The effective thermal conductivity λeff   represents the overall thermal conductivity of solid particles, air in the internal pores of pellets and pellet-to-pellet contact thermal resistance. The commercial pellets were prepared from sawdust pre-dried at temperatures up to 100oC to low moisture contents of 4-8% [119]. The conditions within the bulk are close to equilibrium and the potential for moisture evaporation is minimal under typical storage conditions. Similarly it is assumed that the rate of moisture diffusion can be lumped in an effective thermal conductivity λeff. Table 3.2 shows that the moisture evaporation heat is small compared to the total energy change during experiment. Mohsenin [82] also stated that the effect of moisture migration on estimated thermal conductivity was minimized by experimenting with initially uniform moisture content over  a short time. The heat source method is one of the oldest transient methods and has been widely employed in the engineering field for measuring the thermal conductivity of agriculture materials. Eq. (3.10) was traditionally solved by an analytical method. By assuming that the radius of the vessel is infinite and the radius of the heater is negligible, Eq. (3.10) can be solved [120] to give  2 0 0 exp( ) 2 eff Q x T dx T x        (3.11) where 0 0 02Q r q  is the heating rate per unit length of the heater , w/m; 2 / eff P r t C      with t being the time, s; r0 the radius of the heater, m; q0 the heating rate of the heater source per unit area, W/m2; T0 the initially temperature of the pellets, k; λeff the effective thermal conductivity, W/(m k) Hooper et al. [121] and Ingersoll et al. [120] expanded Eq. (3.11) into a series solution, and obtained an approximate equation using the leading-order terms with the assumption of a negligible product of r and 1/2 1 ( ) 2 t [82, 120-122]. 49  1/20 0 0 0 1 ln( ) 2 2 2 4 eff eff eff Q B Q Q T T r t  (3.12) where B is a constant, α is the thermal diffusivity defined by  λeff/(ρCp). Using this approximate method, a linear relationship between T and ln(t) is obtained with a slope S=Q0/4πλeff. This approximate method is valid only when the measuring point is close to the heater and does not work well for large particles such as pellets. Therefore, instead of the approximate method, we solve Eq. (3.10) numerically using MATLAB directly and apply a least square curve-fitting method to estimate the thermal properties for pellets at several levels of moisture contents. The MATLAB PDE solver, pdepe, is used to solve this parabolic partial differential equation (PDE) Eq. (3.10) in the one space variable r and time with the following boundary conditions B.C.1 0r r , 0 /effq T r   B.C.2    r R , / 0T r   I. C.      0t  , 0T T where r0 is the radius of the center heater, meter; q0 is the heating rate per area of the electric heater, W/m2; T0 is the initial temperature, K; R is the radius of the vessel, m. The pdepe solver converts the PDEs to ordinary differential equations (ODEs) using a second- order accurate spatial discretization based on a fixed set of user-specified nodes. The time integration is done with ode15s, which is especially designed for stiff problems based on the numerical differentiation formulas (NDFs). As the solution, temperatures are expressed as a function of radial coordinate and time for a given initial temperature T0, thermal conductivity λeff and specific heat capacity CP. The thermal properties of wood pellets CP and λeff were determined by minimizing the residue of Eq. (3.13) between the experimental data and the numerical results by a fminsearch function in MATLAB, which uses the Nelder-Mead simplex algorithm as described in Lagarias et al. [123].  150 2 1 Residue ( ) nt C a nt T T       (3.13) 50 where TC is the temperature from the numerical results and Ta is the experimentally measured temperature. nt is the time count (number of measurements); the time interval between each measurement was 1 minute. We took the average of temperature data in the first 600 s as the initial temperature T0 in order to eliminate the fluctuation errors related with the thermocouples. The time history of the temperature data from 0 to 9000 s was used for curve-fitting.  No liquid water was observed after 9000 s heating, so that the moisture condensation heat was neglected. Also, the evaporation heat was neglected since it is small over the range of 9000 s heating (30-50oC). 3.6. Results and discussion Four repeated experiments were carried out to examine the confidence levels and errors of the proposed fitting method using pellet A with 4.6 % moisture content. The results are listed in Table 3.3, which show that the thermal conductivity and heat capacity values estimated by the proposed method are reproducible with low standard deviations.  Figure 3.3 shows that the temperature profiles measured and fitted by the numerical and analytical methods are almost identical within the first 9000 s. even though different outer boundary conditions were used in the numerical (insulated at the wall) and analytical approaches (infinite radius), suggesting that the outer boundary condition is less significant within the first 9000 s of heating). Both continuous models (numerical and analytical) fitted reasonably well to the experimental results. The fluctuation of the experimental data might be caused by the tolerance of the thermocouple and the heterogeneity of the material. Table 3.3 Reproducibility of experimental data-fitted thermal conductivity and heat capacity  Trial 1 Trial 2 Trial 3 Trial 4 Mean  and Variance (at 95% confidence limit)1 λeff, W/(m K) 0.152 0.161 0.153 0.155 0.155±0.008 CP, kJ/(kg K) 1.223 1.141 1.114 1.213 1.173±0.084 Note: 1 Mean and Variance were calculated base on equation x=average(x)+stn-1,α/sqrt(n) [124] 51 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 22 24 26 28 30 32 34 36 38 Time(s) T e m p e ra tu re (o C )   Experimental Results Numerical Results Analytical Results  Figure 3.3 Typical time variation of the temperature profiles from experiments and fitted by  numerical and analytical method. Note: T1, T2 and T3 are the temperature profiles at the tip of thermocouples T1, T2 and T3, respectively (as shown in Figure 2.1). 3.6.1. Effect of the pellet sizes Effect of pellet size on the effective thermal properties was investigated by sieving pellet A (as received) into three size fractions, with the experimental results listed in Table 3.4. The packing porosity, 𝜀, is seen to increase slightly with the mean pellets size for randomly packed bed of samples, and the effective thermal conductivity remains more or less the same. Wood pellets are regularly shaped and highly condensed material, so that the porosity range of the packed bed is relatively narrow and the porosity effect on effective thermal conductivity is not significant in this narrow range. The results in Table 3.4 also show that specific heat capacities are almost independent of the pellet size.   T1 T2 T3 52 Table 3.4 Effect of pellet size on the effective thermal properties (3 repeats) Sizes Bulk density, kg/m3 Porosity,  λeff, W/(m K) CP, kJ/ (kg K) Small (<6.7mm) 670 0.422 0.168±0.010 1.109±0.102 As received 665 0.426 0.160±0.008 1.163±0.080 Large (>6.7mm) 660 0.430 0.164±0.009 1.023±0.103 Note: Small and Large pellets are obtained from as received pellets using a sieving machine. 3.6.2. Moisture content effect on thermal conductivity Experimental results of the thermal properties of pellets with moisture content ranging from 1.4 to 9 % w.b. are given in Tables 3.5 to 3.7. Thermal conductivity of pellet A in Table 3.5 increased from 0.146 to 0.192 W/(m K) when the moisture content increased from 1.4 to 9 % w.b., while the porosity remained close to a constant. The same trend was recorded in Table 3.6 for pellet B over the same range of moisture contents as pellet A. However, it is noted that the porosity of pellet B ranged from 0.418 to 0.463, and it is seen that the lowest effective thermal conductivity of 0.137 was obtained for pellets with a porosity of 0.463, compared to 0.177 for samples with a porosity of 0.418, although the later had a higher moisture content than the former. The results show that the experimental values are larger than the value of 0.0986 W/(m K) reported by Gupta [86] for dried softwood particles (bulk density of 360 kg/m3), and this is partly related to the higher particle density of pellets. Tables 3.5 and 3.6 also show that the effective thermal conductivity values predicted using equations reported in the literature and listed in Table 1.4. The Kollman equations [85] predicts the dry material thermal conductivity well but cannot account for the moisture content effect on the thermal conductivity. In the Kollman equation, the tested wood moisture content ranged from 5% to 35%. However, the maximum moisture content in wood pellets is only 9-10% (wet based) and pellets will break down into sawdust if moisture content exceeds 12%. The predicted values from the Fasina equation [46], which was developed from alfalfa pellets data, are seen to be smaller than our measured wood pellets values although the predicted trend of the moisture effect is similar.     53 Table 3.5 Thermal conductivity of pellet A obtained from experiments and predicted from model equations M.C., % w.b. Bulk density1, kg/m3 Porosity, λeff, W/(m K) Experiments λeff, W/(m K) Kollman2 λeff, W/(m K) Fasina3 λS, W/(m K) Parallel4 1.7 650 0.421 0.146±0.014 0.154 0.063 0.233 4.6 665 0.426 0.155±0.006 0.154 0.087 0.250 6.5 665 0.426 0.160±0.008 0.154 0.102 0.259 8.5 670 0.433 0.173±0.010 0.154 0.119 0.284 9.0 675 0.434 0.192±0.002 0.154 0.123 0.319 Note: 1 Bulk density in this table is wet based. 2 Kollman [85] effective thermal conductivity is calculated using equations  0.0002 0.024dry dry    and 2(1 0.00015 )eff dry M   .  3 Fasina [46] effective thermal conductivity is calculated using equation 0.049 0.0082eff M   . 4 λS parallel is the single pellet thermal conductivity calculated using Eq.(3.7).  Table 3.6 Thermal conductivity of pellet B obtained from experiments and predicted from literature equations M.C., % w.b. Bulk density1, kg/m3 Porosity, λeff,W/(m K) Experiments λeff,W/(m K) Kollman2 λeff,W/(m K) Fasina3 λS, W/(m K) Parallel4 1.4b 590 0.463 0.139±0.008 0.140 0.060 0.235 4.6a 690 0.418 0.177±0.010 0.156 0.087 0.285 6.5a 690 0.430 0.181±0.012 0.153 0.102 0.294 7.3a 700 0.427 0.187±0.014 0.154 0.109 0.306 8.8b 650 0.443 0.185±0.018 0.143 0.121 0.311 Note: a is conditioned (or as received) from pellet B1. b is conditioned (or as received) from pellet B2. 1 Bulk density in this table is wet based. 2 Kollman [85] effective thermal conductivity is calculated using equations 0.0002 0.024dry dry    and 2(1 0.00015 )eff dry M   .  3 Fasina [46] effective thermal conductivity is calculated using equation 0.049 0.0082eff M   . 4 λS parallel is the single pellet thermal conductivity calculated using Eq. (3.7).       54 The single pellet thermal conductivity for each sample is also estimated following the parallel distribution model, which accounts for the effect of packing porosity, and listed in Table 3.5 and Table 3.6. The single pellet thermal conductivity is seen to range from 0.233 to 0.319 W/(m K) and increases with increasing the moisture contents of the pellets. Figure 6 shows a linear regression of the moisture content and the single pellet thermal conductivity using the experimental data from both pellet A and pellet B. 0.219 0.01s M                                                     (3.14) where λs is the single pellet thermal conductivity, in W/(m K); M is the moisture content, in percentage (wet basis). At the tested temperature range (25oC-50oC), the moisture evaporation is limited according to Appendix A so that the moisture content effect on the effective thermal conductivity can be estimated following the solid-liquid layer system models described in Eq. (3.7) and Eq. (3.8). From Figure 3.4, we can see that the moisture content affects the single pellet thermal conductivity much more than predicted by both the parallel and serial distribution models. Combining the parallel model to account for the porosity effect (for calculating single pellet thermal conductivity) and Eq. (3.14), Eq. (3.15) is obtained by least-square regression for predicting the effective thermal conductivity of packed wood pellets.  (0.219 0.01 ) (1 ) 0.027eff M         (3.15)  where λeff  are the thermal conductivity of the packed bed and the single pellets, W/(m K) and ϕ is the porosity in fraction; M is the moisture content in percentage, wet based. 55 0 2 4 6 8 10 0.15 0.20 0.25 0.30 0.35 0.40    Pellet A  Pellet B  Eq. (3.14)  Parallel distribution model  Serial distribution model S in g le  p e lle t th e rm a l c o n d u c ti v it y  ( W /( m k )) Moisture content (%)  Figure 3.4 Comparison of experimental results (pellet A and pellet B) and predicted single pellet thermal conductivity using Eq.(3.14), parallel distribution model and serial distribution model. Note: Parallel distribution model was calculated using equation , / 100 (1 / 100) s water s dry M M      Serial distribution model was calculated using equation , 1 100 1 100 s water s dry M M       where 0.626 water    W/(m K) and ,s dry   was determined to be 0.255 and 0.259 W/(m K) for parallel distribution and serial distribution respectively.       56 3.6.3. Moisture content effect on specific heat capacity Table 3.7 shows the experimental results of specific heat capacities for both pellet A and pellet B. The specific heat capacity of wood pellets ranges from 1.074 to 1.253 kJ/(kg K) in the test range. This result is similar to those reported by Adl-Zarrabi, 1.07 kJ/(kg K) for dry wood and 1.38 kJ/(kg K) for M.C 9.5 % wood [83]. When the average specific heat capacity (Cp,M) of water (37 oC) is set to 4.178 kJ/(kg K) [112], the bulk pellets with moisture can be treated as a mixture of dry pellets and water, so that the specific heat capacity becomes , (1 /100) 0.042 P P dry C M C M   . (3.16) Using the experimental results of both pellet A and pellet B, the specific heat capacity of dry wood pellets (CP,dry) was found to be 1.01±0.05 kJ/(kg K) and Eq. (3.16) became 1.01 0.032 P C M  . (3.17) The specific heat capacity of dry wood pellets, 1.01±0.05 kJ/(kg K), is also close to those predicted from Gupta [86], Kollman [85] and Fasina [46] equations, which are 1.17, 1.31 and 1.08 kJ/(kg K) for wood particle, general wood and alfalfa pellets, respectively.  Table 3.7 Specific heat capacities of wood pellets obtained by for pellet A and pellet B Pellet A Pellet B M.C., % w.b. CP, kJ/(kg K) M.C., % w.b. CP, kJ/(kg K) 1.7 1.099±0.150 1.4 1.073±0.052 4.6 1.173±0.092 4.6 1.012±0.090 6.5 1.163±0.080 6.5 1.190±0.006 8.5 1.023±0.103 7.3 1.219±0.090 9.0 1.185±0.140 8.8 1.253±0.080  57  3.7.  Conclusions The effective thermal conductivity and specific heat capacity of wood pellets were experimentally determined using a modified line heat source experimental method [46, 121] and a Matlab curve- fitting data analysis method in this Chapter. Two types of wood pellets produced in British Columbia were investigated and the following conclusions are obtained. 1. The numerical curve-fitting method which was developed to determine the thermal properties, including effective thermal conductivity and specific heat capacity of wood pellet provides reproducible results with small standard deviations. 2. Moisture content has a significant effect on effective thermal conductivity and specific heat capacity in the range of 1.4 to 9% w.b. studied. Effective thermal conductivity is much less sensitive to bulk density than to moisture content and insensitive to the size of the pellets. 3.  Empirical correlations were developed based on the experimental results to predict the effective thermal conductivity and specific heat capacity of pellets as a function of  their moisture content and bulk density. 58  4. Measurements of self-heating using isothermal calorimetry1 Self-heating of bulk pellets at low temperatures is a slow process that can only been observed in large vessels. The heat release is too small to be detected using normal laboratory equipment. This Chapter presents the results of   heat generation rates at low temperatures (30-50 oC)  measured using an isothermal calorimetry method , and the kinetic parameters, including activation energy E, reaction enthalpy ∆rh, and pre-exponential factor A, obtained from regression analysis of experimentally measured heat production rate data. The factors affecting the self-heating reactions, including the temperature, humidity, oxygen supply in the storage as well as the moisture contents, type and age of the pellets are discussed. . 4.1. Materials Softwood pellets produced in British Columbia (BC) were used in this research. Pellets with designation A, B and C were manufactured from 3 different wood pellet plants in BC. Pellet A and pellet B were made of whitewood, mostly pine; pellet C was made of a mixture of pine and a small amount of bark.  Pellet A and pellet C were manufactured by Princeton Co-Generation Crop. (Princeton, BC, Canada). Pellet B was manufactured by Premium Pellet Ltd. (Vanderhoof, British Columbia Canada). The moisture content was measured following the ASABE standard S358 using the precision oven model 70D (Thermo Fisher Scientific Inc., 308 Ridgefield Court Asheville, NC 28806). Table 4.1, 4.2 and 4.3 list the properties of the tested 3 groups of pellets. Group 1 was composed of as received pellets and some aged pellets. In this chapter, age of the pellets is defined as a function of the time and temperature of the storage. For instance, sample 5 was stored at ambient temperature for a year after it was produced thus its assumed age is 1 year old. Sample 4 was conditioned from sample 5 using the ESPEC humidity chamber model LHU-113 (ESPEC North America, Inc., 4141 Central Parkway Hudsonville, MI 49426). Sample 9 was dried using a  1 A version of this chapter was submitted for publication. Guo, W., Trischuk, K., Bi, X., Lim, CJ., Sokhansanj, S.. Measurements of wood pellet self-heating kinetic parameters using isothermal calorimetry. 59 freeze dryer to compare with the oven dried pellets for the purpose of studying the effect of treatment temperature. Table 4.1 Properties of group 1 samples (as received pellets and aged pellets) Sample 1 2 3 4 5 6 7 8 9 pellets B B B A A A B C B M.C. (w.b.) a 4 % 8 % 0 % 9 % 7 % 0 % 6 % 4 % 0 % Cp , kJ/(kg.K) 1.2 1.3 1.1 1.3 1.3 1.1 1.2 1.2 1.1 Age b & condition Fresh Fresh Oven dried 1 year 1 year Oven dried 2 months 2 months Freeze dried a   M.C. stands for moisture contents of the pellets, wet based. b Ages were defined as a function of the time and condition of storage, assuming that pellets age will increase  when they are in storage, unless they are stored in cold room (temperature under 0oC).  For instance, sample 5 was stored at room temperature for a year after received from the factory thus its assumed age is 1 year old.  Table 4.2 Properties of group 2 samples (stored at room temperature 23oC) Sample R1 R2 R3 R4 R5 R6 R7 Pellets B B B B B B B M.C. (w.b.) a 4 % 4 % 4 % 4 % 4 % 4 % 4 % Conditioned time 38 days 28 days 21 days 17 days 11 days 6 days 3 day a   M.C. stands for moisture contents of the pellets, wet based.  Table 4.3 Properties of group 3 samples (stored under 50oC) Sample F1 F2 F3 F4 F5 F6 Pellets B B B B B B M.C.  (w.b.) a 4 % 4 % 4 % 4 % 4 % 4 % Conditioned time  6 days 5 days 4 days 3 days 2 days 1 day a   M.C. stands for moisture contents of the pellets, wet based.  In order to study the age effect, Group 2 (R1-R7) were fresh pellets conditioned at room temperature (23oC) and Groups 3 (F1-F6) were fresh pellets conditioned in a oven (50oC) for different length of time. To prepare the Group 2 samples, 7 fresh pellet samples, 50 g each, were taken from the refrigerator and put into the conditioning cabinet at room temperature.  R7 sample aged for 3 days was  the sample  taken out from the cabinet after 3 days, put in a sealed bag, and stored  in the refrigerator until the test.. R6 (aged for 6 days) was samples removed from the cabinet 60 after 6 days , put in a sealed bag,  and put back into the refrigerator. The sample removed from the cabinet after 38 days and stored in refrigeration is marked as R1.  The preparation of Group 3 followed the same procedure. 4.2.  Self-heating kinetics and theories According to Kuang et al. [33-35], the decomposition reaction of wood pellets can be assumed to be first order reactions with main products as CO2, CO and CH4 and heat (q).  The overall reaction can thus be written as Pellets + O2 → CO2 + CO + CH4 + q With the assumption that the pellets decomposition reactions can be separated into oxygen dependent reactions and oxygen independent reactions, the heat (q) from two groups of reactions can be expressed as  ox noxq q q   (4.1) where qox is the heat released from oxygen dependent reactions, in J; qnox is the heat released from oxygen independent reactions, in J. 4.2.1. Oxygen independent reactions The test is conducted in a closed ampoule in a constant temperature and starting with atmospheric air. As the oxygen is consumed by the reactions, the heat released from oxygen-dependent reactions decreases while the heat released from oxygen independent reactions remains the same. At the late stage of the tests, there is no or little oxygen left in the ampoule with most of the oxygen consumed by the reactions. The measured heat can be attributed to the heat release from oxygen independent reactions. The heat released from the oxygen independent reactions (qnox) can thus be obtained, as shown in Figure 4.1. 61 4.2.2. Oxygen dependent reactions Assuming oxygen dependent reactions can be lumped as a single first order reaction with respect to oxygen and there is abundant amount of biomass, the reaction rate of oxygen dependent reactions is given by  2 2 O A O dC r kC dt     (4.2) where rA is the reaction rate base on O2 at time t, mol/( m 3.s); k is the reaction rate constant, in s-1; CO2 is the concentration of O2 at time t, in mol/m 3.  0 0(1 ) dX C kC X dt    (4.3) where X is the fraction of conversion of the oxygen dependent reaction, C0 is the initial concentration of O2 , in mol/m 3. Heat release from oxygen dependent reactions can be expressed by  0ox rq C VX H   (4.4) where V is the volume of the container, in m3; ∆rH is the reaction enthalpy (reaction heat) of the oxygen dependent reactions, in J/mol. For the heat release rate, inserting Eq. (4.4) for X into Eq. (4.3), one obtains  0 0( ) ( ) ox ox r r ox r dq q k H C V k C V H q dt H         (4.5) Solving Eq. (4.5), we get  1 ktoxdq Ae dt   (4.6) 62 where A1 is a heat release rate coefficient, in J/s ; k is the reaction rate constant, in s -1. From the TAM isothermal calorimeter, we can get dqox/dt as a function of time (t). The difference between dqox/dt value of experiments and Eq. (4.6) using different A1 and k can be calculated.  exp .(06) exp ( ) / ( )ox ox ox Eq dq dq dq Error dt dt dt    (4.7) By minimizing the errors, a best fitted values of A1 and k can be obtained. Reaction rate constant (k) is normally a strong function of temperature. Arrhenius equation has been commonly used to describe the temperature dependence of the reaction rate constant (k).  / 0( ) E RTk T A e  (4.8) where A0 is the pre-exponential factor or frequency factor; E is the activation energy, in J/mol; R is gas constant, 8.314 J/(mol∙K); and T is the absolute temperature, K. Assuming that the reaction enthalpy remains constant over the tested temperature range and all the oxygen is consumed in the reaction (when time is infinite), the reaction enthalpy of oxygen dependent reactions can be calculated using Eq. (4.4) and Eq. (4.6)  1 0 s r A M H kVC    (4.9) where V is the volume of the container, in m3; ∆rH is the reaction enthalpy (reaction heat) of the oxygen dependent reactions, in J/mol; A1 is a heat release rate coefficient, in J/s. Ms is the sample weight during the test. 4.3. Results of the heat release rate The heat release rate of group 1, 2 and 3 were tested using TAM air isothermal calorimeter at 30, 40 and 50oC.  For each test, the heat flux from the samples (wood pellets) was measured at certain constant temperature at which the TAM calorimeter was previously calibrated. The values obtained 63 from heat release rate measurements were divided by the sample weight in order to get the heat release rate per unit weight. 4.3.1. Temperature effect on the heat release rate Figure 4.1 shows a typical plot of heat release rate for sample 2 from group 1 (fresh pellet B with 8% w.b. moisture content) in a closed ampoule using the TAM air isothermal calorimeter at 30, 40 and 50oC. The heat release data indicate that the test temperature has a significant impact on the heat release rate. The heat release rate is higher at the elevated temperature. A maximum self-heating release rate of 0.33 mW/g was observed at 50oC. Initially, the heat release rate reaches a peak value (0.33 mW/g at 50oC and 0.04 mW/g at 30oC) and then decreases after the peak period. We speculate that the heat release rate decreased with time after the peak was reached due to the decrease in the oxygen concentration, and finally reached a constant value after 50 hours (1.8x105 s), which corresponds to the heat release rate of the oxygen independent reactions. All tests were measured in two repeats. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Time (s) H e a t re le a s e  r a te  ( m W /g )   50oC 40oC 30oC  Figure 4.1 Measurement results of sample 2 (Fresh pellet sample B with 8% w.b moisture content) in closed ampoule using a TAM air isothermal calorimeter at 30, 40 and 50oC Oxygen independent reactions dominate Oxygen dependent reactions dominate 64 At the initial period, where the heat production from the pellets appears to be the highest, the heat release must have been dominated by oxygen dependent reactions; at the end, where the heat production decreases and the oxygen is consumed, the heat is mainly released from oxygen independent reactions. This observation is similar to Kuang 's[34] research finding in the off-gassing emission research. They observed a high reaction rate at the initial period, slowly decreased reaction rate after the initial period and an equilibrium stage at the end where the CO2 and CO concentrations stopped increasing. Kuang’s research also showed that oxygen levels showed an exponential decrease, similar to the decrease in the heat release rate observed from TAM experiments. 4.3.2. Moisture content effect on heat release rate Sample 1 and sample 2 (group1) are freshly made to study the moisture effect on self-heating. The pellets were made by the same manufacturer, using the same material and same procedure except the pellets had different moisture content (sample 1 had 4% moisture content and sample 2 had 8% moisture content). Sample 3 has oven dried from sample 1 and sample 9 was freeze dried from sample 1. Both oven dried pellets and freeze dried pellets had 0% moisture content. Because the 50oC test showed the highest self-heating rate in the tested range, the results of 50oC measurement were used to study the effect of moisture content on the heat release rate. 65 0 0.5 1 1.5 2 2.5 3 x 10 5 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Time (s) H e a t re le a s e  r a te  ( m W /g )   Sample 1(B,fresh, M.C 4%) Sample 2(B,fresh, M.C 8%) Sample 3(B,oven dried, M.C 0%) Sample 9(B,fresh,freezing dried,0%)  Figure 4.2 Measurement results of pellets B with different moisture content (sample 1, 2, 3, 9)  Figure 4.2 shows measurement results at 50oC for sample 1, 2, 3 and 9 (group1). Sample 1 (4% moisture content) and sample 2 (8% moisture content) gives almost the same amount of heat release rate. Almost no heat is detected in sample 3 (oven dried 0% moisture content) ampoule, while the freeze dried pellets (sample 9) give the same amount of heat as the fresh pellets (sample 1 and 2). Therefore, it appears that the moisture content of pellets within these low ranges does not have significant effects on the self-heating rate whereas storage temperature has an important impact on the self-heating rate. 4.3.3. Age effect on the heat release rate The age of the wood pellet is defined as the number of storage days before a performing self-heating test. The age of fresh pellets after its manufacture is defined as 0. In this section, sample 1, sample 3, sample 5, sample 7 and sample 8 were used to illustrate the age effect on the heat release rate and the tests results of 50oC were shown in Figure 4.3. All the ages indicated in Figure 4.3 are ages under room temperature. Fresh pellets (sample1) give the highest heat release rate; two-month old pellets 66 (samples 7 and 8) have lower heat release rate and oven dried pellets (sample 3) or pellets with an age of more than one year (sample 5) emit almost no heat. Therefore, age plays an important role during the self-heating process of wood pellets storage. The self-heating rate pellets from different manufactures (samples 4, 5, 6 and 8) also were investigated. No significant impact of pellets type was observed on the self-heating rate of wood pellets. For instance, Figure 4.3 illustrates the heat release rate of pellet A (sample 5) and pellet C (sample 8), and the test results are close to that test results of pellets B. Therefore, in this research, the kinetics parameter determination was mainly focused on the age effect and temperature effect. The group 2 and group 3, which were conditioned at 23oC and 50oC respectively for different days, were tested at 30, 40 and 50oC to further study the effect of age and temperature effect and to determine the self-heating kinetics. 0 0.5 1 1.5 2 2.5 3 x 10 5 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Time (s) H e a t re le a s e  r a te  ( m W /g )   Sample 1(B,fresh, M.C 4%) Sample 7(B,2 month old, M.C 6%) Sample 8(C,2 month old, M.C 4%) Sample 3(B,oven dried,0%) Sample 5(A,one year old,7%  Figure 4.3 Measurement results for group 1 samples at 50oC in closed ampoules using a TAM air isothermal calorimeter. ‘A’, ‘B’ and ‘C’ stand for pellet type A, B and C, respectively.  67 4.4. Heat from oxygen independent reactions The heat released from oxygen independent reactions can be calculated from the last period of tests where the oxygen dependent reactions are almost complete. The results show that oxygen independent reactions are not sensitive to age, moisture contents or pellets types. Heat release from oxygen independent reactions is rarely temperature dependent in our study range. To reduce the system noise and errors, all test results from group 1, group 2 and group 3 were averaged to obtain an average heat release rate for the oxygen independent reactions over the test range, with the resulting values of 0.0026± 0.0009 mW/g obtained. 4.5. Kinetic and heat release rate parameters of oxygen-dependent reactions Pellets samples from group 1, group 2 and group 3 are used in this section to study the kinetic parameters of oxygen dependent reactions. Experiments were conducted at temperature 30oC, 40oC and 50oC and the heat release rate from each sample was obtained at each temperatures. Assuming that both the reaction rate constant k and coefficient A1 remained constants over the range of tests, a curve-fitting method was developed to estimate the value of k and A1 by solving Eq. (4.6) using Matlab, and by minimizing the error between experimental data and calculated data using  Eq. (4.7). The estimated values for k and A1 will be discussed in the following section. Moreover, the activation energy of the fresh pellets and the reaction enthalpy of the oxygen dependent reaction will also be discussed. 4.5.1. Reaction rate constant k As previously discovered in Figure 4.2, the fresh pellets give out the highest heat release rate regardless the pellets moisture content. In this section, the reaction rate constant k and the activation energy of fresh samples (samples 1, 2 and 9) are calculated and listed in Table 4.4. All tests were measured in two repeats. Similar to the observation in Figure 4.2, reaction rate constant k appears to be strongly dependent on temperature where a higher k values are found at a higher temperature. For fresh pellets (age zero), the average k value is estimated to be 4.77*10-6 s-1, 8.10*10-6 s-1 and 1.91*10-5 s-1 for 30oC, 40oC and 50oC, respectively.  68 Table 4.4 Reaction rate constant and activation energy of fresh samples with different moisture contents Temperature Sample 1 (Pellets B, M.C. 4%) Sample 1 (Pellets B, M.C. 4%) Sample 2 (Pellets B, M.C. 8%) Sample 2 (Pellets B, M.C. 8%) Sample 9 (Pellets B, M.C. 0%) Sample 9 (Pellets B, M.C. 0%) k (30 oC) *105, s-1 0.489±0.003 0.549±0.001 0.415±0.002 0.543±0.002 0.493±0.002 0.372±0.010 k (40 oC) *105, s-1 0.941±0.003 0.723±0.001 0.748±0.001 0.987±0.003 0.770±0.002 0.688±0.005 k (50 oC) *105, s-1 2.094±0.013 2.035±0.012 2.176±0.009 1.795±0.012 1.704±0.012 1.676±0.010 E (kJ/mol) 56.0±14.3 57.9±14.1 55.7±12.1   4.5.1.1. Activation energy and age The activation energy was calculated using Arrhenius equation (Eq.(4.8)), with the results also listed in Table 4.4. The activation energy of pellets ranges from 50 to 61 kJ/mol, slightly lower than those reported in the literature for woody materials, 60-80 kJ/mol [24, 26, 92]. This lower activation energy could be caused by the bacteria present in the fresh pellets, which act like a catalyst to lower the activation energy and promote the decomposition process [23, 41]. -1 0 1 2 3 4 5 6 7 8 50 55 60 65 70 75 80 85 90   A c ti v a ti o n  e n e rg y  ( k J /m o l) x 1   (days)  Figure 4.4 Effect of aging on activation energy for group 2 pellets (stored at 23oC).  x1 is the days stored at 23oC before the test. 69 From the previous discussion, we know that pellet's age is a significant factor affecting the reaction rate. In order to further study the age effect and storage temperature effect on the reaction rate constant and activation energy, the self-heating of group 2 and group 3 samples were tested. They were conditioned at 23oC and 50oC for different times, respectively. Figure 4.4 shows the effect of aging on the activation energy using group 2 pellets (conditioned at 23oC). The activation energy increases with age (conditioned days). It is seen that the activation energy does not change much at the first 4 days of conditioning but increases rapidly after conditioned for more than 4 days and reaches 90 kJ/mol at the 7 days of conditioning. 4.5.1.2. Correlation of age and temperature The reaction rate constants k of group 2 and group 3 were calculated in determine the effects of pellet aging condition and test temperature. The reaction rate constant k at test temperatures of 30oC, 40oC and 50oC for group 2 samples (conditioned at 23oC) are listed in Figure 4.5. At 50oC, the reaction rate constant is the highest and the lowest is at 30oC, which is identical to what is observed in group 1 pellets. As the aging increases, the reaction rate constant decreases rapidly at a test temperature of 50oC but decreases much slower at test temperatures of 40oC and 30oC. In another word, self-heating reaction releases heat faster at higher temperatures and the pellet ages faster at higher temperatures as well. During the test of group 3 pellets (conditioned at 50oC), the heat release was detected at 50oC test, but  no heat release was detected in 30oC and 40oC tests. The reaction rate constant k (50oC), the k value at 50oC, of group 2 and group 3 pellets were listed in Figure 4.6. The reaction rate constant k (50oC) for group 2 pellets (conditioned at 23oC) fits well to an exponential relationship with the pellet age, with a coefficient of determination R2 = 0.934, as shown in Figure 4.6.  5 1(50 ) 2.063exp( 0.014 ) 10 ok C x     (4.10) where k(50oC) is the reaction rate constant at 50oC.  x1 is the number of days pellets were stored  at 23oC before the test. 70 -5 0 5 10 15 20 25 30 35 40 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2    k(30 o C)  k(40 o C)  k(50 o C) R e a c ti o n  r a te  c o n s ta n t *1 0 5  ( s -1 ) x 1  (days) Eq. (12)  Figure 4.5 Comparison of experimental reaction rate constant data k(T) and Eq. (12) for group 2 pellets (stored at 23oC). x1 is the days stored at 23oC before the test. Note: T is the temperature of test. The k(50oC) values of group 3 pellets (stored at 50oC) decrease much more expeditiously with aging (storage days) compared to group 2 pellets (stored  at 23oC). The reaction rate constant, k(50oC), for group 3 pellets (stored at 50oC) also fits well to an exponential relationship with a coefficient of determination (R2) of 0.915, as shown in Figure 4.6, given by Eq. (4.11),.  5 2(50 ) 1.894exp( 0.255 ) 10 ok C x     (4.11) where k(50oC) is the reaction rate constant at 50oC, and  x2 is the days conditioned at 50 oC before the test.  71 -5 0 5 10 15 20 25 30 35 40 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2    Group 2  Group 3  Eq. (10)  Eq. (11) R e a c ti o n  r a te  c o n s ta n t *1 0 5  ( s -1 ) x 1 (days)  Figure 4.6 Comparison of experimental reaction rate constant data and Eq. (4.10) and Eq.(4.11) for 50oC tests of group 2 pellets (stored at 23oC) and group 3 (stored at 50oC), respectively. Note: x1 is the days stored before the test.  Comparing Eq. (4.10) with Eq.(4.11), the reaction rate constants of group 3 pellets (stored  at 50oC) decrease about 19 times faster than group 2 pellets (stored at 23oC). This relation can be used to combine Eq. (4.10) and Eq.(4.11) into one equation in order to make a unified prediction of the self- heating rate. Combining Eq. (4.8) and Eq. (4.10), with an average activation energy of 56.55kJ/mol for fresh pellets, k(T) for group 2 pellets (conditioned at 23oC) can be predicted by  4 1 1 6801.78 ( , ) 2.883exp( 0.014 ) 10k x T x T      (4.12) where k(x1,T) is the reaction rate constant for tests at temperature T for pellets stored in 23 oC for x1 days before the test. Figure 4.5 show that Eq. (4.12) fits the experimental data well. 72 4.5.2. Heat release coefficient A1 Heat release coefficient A1 was determined along with the reaction rate constant k for pellets of group 1, group 2 and group 3. Table 4.5 lists the heat release rate coefficient A1 for group 1 pellets of different moisture contents. Similar to the reaction rate constant, the heat release rate coefficient A1 strongly depends on the test temperature and is the largest at 50 oC and smallest at 30oC. Even though it is also slightly sensitive to the moisture content of pellets, the difference before the highest number (8% M.C) and the smallest number (0% M.C) is really small. Averaged over six samples at each temperature, A1 values are obtained as 0.028, 0.065 and 0.247 mW/g for 30, 40 and 50 oC, respectively. Table 4.5 Heat release coefficient A1 data (mW/g) of Group 1 fresh pellet samples with different moisture contents. Temperature Sample 1 (Pellets B, M.C. 4%) Sample 1 (Pellets B, M.C. 4%) Sample 2 (Pellets B, M.C. 8%) Sample 2 (Pellets B, M.C. 8%) Sample 9 (Pellets B, M.C. 0%) Sample 9 (Pellets B, M.C. 0%) A1(30oC)*102 2.415±0.011 2.308±0.003 2.670±0.008 2.605±0.008 1.558±0.003 2.157±0.002 A1(40oC)*102 10.197±0.022 4.108±0.004 7.847±0.015 4.804±0.005 6.984±0.011 5.023±0.025 A1(50oC)*102 22.844±0.098 25.534±0.107 29.095±0.087 25.593±0.118 24.636±0.118 20.792±0.087  Same as reaction rate constant equations, the heat release coefficient A1 (50 oC), A1 value at 50 oC, fits well to exponential expressions given by Eq. (4.13) and Eq.(4.14), for group 2 and group 3 pellets, respectively. Figure 4.7 compares the experimental A1 values with Eq. (4.13) and Eq. (4.14) for group 2 pellets (conditioned at 23oC) and group 3 pellets (conditioned at 50oC). Group 2 1 1(50 ) 0.276exp( 0.0256 ) oA C x        (R2=0.969) (4.13) Group 3 1 2(50 ) 0.254exp( 0.640 ) oA C x          (R2=0.982) (4.14) where A1(50 oC) is the heat release coefficient at 50oC, x1 is the days conditioned at 23 oC before the test and x2 is the days conditioned at 50 oC before the test. 73 -5 0 5 10 15 20 25 30 35 40 0.00 0.05 0.10 0.15 0.20 0.25 0.30    Group 2  Group 3  Eq. (14)  Eq. (15) A 1  ( m W /g ) x 1  (days)  Figure 4.7 Comparison of heat release coefficient A1 data and Eq. (13) and Eq. (14) for tests at 50oC using group 2 pellets (stored at 23oC) and group 3 pellets (stored at 50oC) respectively. x1 is the days stored before the test. Assuming that A1 follows the same trend as reaction rate constant k, the following expression is found to fit well of all data obtained from different test temperatures for pellets conditioned at 23oC.  15 1 1 1 11637.80 ( , ) 1.151exp( 0.0256 ) 10A x T x T      (4.15) where A1(x1, T) is the reaction rate coefficient during test at temperature T for pellets stored at 23 oC for x1 days before the test. As shown in Figure 4.8, Eq. (4.15) fits well with all the experiments data. 74 -5 0 5 10 15 20 25 30 35 40 0.00 0.05 0.10 0.15 0.20 0.25 0.30   A 1 (30 o C) A 1 (40 o C) A 1 (50 o C) A 1  ( m W /g ) x 1  (days)  Figure 4.8 Comparison of heat release coefficient A1 data A1(T) and Eq. (15) for group 2 pellets (stored at 23oC). Note: x1 is the days stored at 23oC before the test. T is the temperature of tests.  4.5.3. Reaction enthalpy Assuming that all oxygen was consumed by reactions at the end of the test, the reaction enthalpy in kJ/mol(O2) can be calculated using Eq.(4.9). To convert the measured heat release rate in mW to the reaction enthalpy in J/mol(O2), the volume of oxygen in the ampoule at the start of tests is needed. For a volume of the closed ampoule of 23.5 ml with 4.5 ml being occupied by 5 grams pellet sample (estimated based on the single pellet density of 1.17 g/ml from a previous study [125], the air volume in the closed ampoule is 19 ml at an initial ambient pressure (with 21% of oxygen), which gives an initial oxygen content of 164 µmol in the ampoule. The heat release rate in mW/mol(O2) is then estimated. The reaction enthalpy is calculated using Eq.(4.9) and experimental results, and the average values of 170, 242 and 464 kJ/mol(O2) at 30, 40 and 50 oC, respectively, were obtained for fresh pellets. The reaction enthalpy is not sensitive to moisture contents of pellets samples. However, they decrease with aging, likely caused by the incomplete consumption of oxygen in the ampoule 75 because the activation energy increases with aging time and the reactions become more difficult to proceed. 4.6. Comparison of self-heating kinetic model with experimental data A comprehensive self-heating kinetic model was obtained by combining the oxygen dependent and independent reactions. As given by Eq. (4.16), the first term on the left represents the heat generation from oxygen dependent reactions and the second term represents the heat generation from oxygen independent reactions. The total heat generation rate is expressed as a function of pellets’ initial age and the storage temperature. Therefore, this model can be used to simulate the self-heating process at different storage temperatures for pellets with different initial ages.  15 1 1 4 1 11637.80 ( , ) 1.151exp( 0.0256 ) 10 6801.78 exp[ 2.883exp( 0.014 ) 10 ] 0.0026 dq x T x dt T x t T             (4.16) where dq/dt is the heat release rate of pellets with initial age x1 at temperature T, in mW/g. T is the absolute temperature in the storage vessel in K, x1 is the initial age in days stored at 23 oC before loaded into the container, and t is the time of test (or reaction time) in s at temperature T. Figure 4.9, Figure 4.10 and Figure 4.11 show the comparison of measured heat release rate data with Eq.(4.16) for tests at 50oC, 40oC and 30oC, respectively.  The results show good agreement between the experimental heat release rate data and those predicted by Eq.(4.16), especially for pellets with a short storage time. 76 0 0.5 1 1.5 2 2.5 3 x 10 5 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Time (s) H e a t re le a s e  r a te  ( m W /g )   Eq. (16) experimental data Fresh Pellets 10days 38 days 2 months 6 months 12 months 12 months  Figure 4.9 Comparison of tested heat release rate data and predicted by Eq. Eq.(4.16) at reaction temperature 50oC for pellets conditioned at 23oC (ambient temperature). The age of pellets at the time of test is shown on the graph. 0 0.5 1 1.5 2 2.5 3 x 10 5 -0.1 -0.05 0 0.05 0.1 0.15 Time (s) H e a t re le a s e  r a te  ( m W /g )   Eq.(16) experimental data Fresh Pellets 17 days 38 days 2 months one years old   Figure 4.10 Comparison of tested heat release rate data and predicted by Eq.(4.16) at reaction temperature 40 oC for pellets conditioned at 23oC (ambient temperature). The age of pellets at the time of test is shown on the graph. 77 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x 10 5 -0.1 -0.05 0 0.05 0.1 0.15 Time (s) H e a t re le a s e  r a te  ( m W /g )   Eq. (16) experimental data Fresh Pellets 2 months one years old  Figure 4.11 Comparison of tested heat release rate data and predicted by Eq. Eq.(4.16) at reaction temperature 30oC for pellets conditioned at 23oC (ambient temperature). The age of pellets at the time of test is shown on the graph. 78  4.7.  Conclusions The self-heating from decomposition of stored wood pellets is experimentally measured using a TAM air isothermal calorimeter, with the kinetic and heat release rate parameters being determined using curve-fitting method in this Chapter. The investigation is based on the test of wood pellets produced in British Columbia, Canada. The following conclusions are reached. 1.Test temperature and pellets age are two major factors impacting the self-heating of tested wood pellets samples. Higher heat release rates are associated with higher test temperatures within the tested range, while older pellets released less heat compare to freshly made pellets. However, the heat release rate is less sensitive to the moisture content and the manufacturer where the pellets are made in this investigation.  2.The activation energy of fresh wood pellets determined in this study ranges from 50-60 kJ/mol. The activation energy also increases with the age of the pellets. The slightly lower activation energy than those published for wood (60-70 kJ/mol ) is attributed to a degree of microbial activities.  3.A kinetic model equation for self-heating heat release rate is developed based the experimental data, which can be used to simulate the self-heating at different storage temperatures for wood pellets at different initial ages.   79  5. Measurements of wood pellets off-gassing in storage1 Off-gassing, occurred simultaneously with the self-heating of wood pellets during storage, is another serious hazard caused by the decomposition of wood pellets. .   Kuang et al. [33-35, 126, 127] found that the off-gases produced from bulk wood pellets are mainly consisted of CO2, CO and CH4, and the effects of storage temperature, humidity and head space were investigated.  In Chapter 4, we found that pellets age played an important role in the self-heating process. This chapter describes the experiments carried out to study the effects of the pellet ages and  the storage container sizes on off-gassing. 5.1. Materials Softwood pellets used in this research were manufactured by Premium Pellet Ltd. (Vanderhoof, BC, Canada). Pellets are in average 6 mm in diameter and 10 to 25 mm in length, which were made from whitewood residues with no bark. Pellets were packaged in bags (40 lbs/bags) and shipped directly to our lab immediately after being produced. They were kept in a cold room (4oC) before use. Experiments were carried out in three sizes of reactors (45, 10 and 2 liters) at 22oC using fresh pellets to examine the effect of reactor sizes and head spaces on wood pellets off-gassing.  The reactors are shown in Figure 2.2, and samples and experimental conditions are listed in Table 5.1a. A group of pellets aged for about 3 months at different storage conditions were used to investigate the age effect on the off-gassing emissions in a 10L reactor. The tests were conducted at 22oC, 40oC and 50oC and the test conditions are listed in Table 5.1b.    1 A version of the chapter will be submitted for publication. Guo, W., Lim, CJ., Bi, X., Sokhansanj, S., Melin, ,S.. Effect of pellets age on off-gassing emission from wood pellets in storage. 80 Table 5.1 Sample information and experimental conditions for off-gassing study a. Sample groups for reactor size and head space effects (fresh pellets tested at 22oC) Sample A1 B1 C1 C2 C13 C14 Reactor Size 45L 45L 10L 10L 2L 2L Head space 0% 25% 25% 0% 0% 25%  b. Sample group C for age and head space effects (aged pellets prepared for tests at 22oC) Sample  C1 C2 C3 C4 C5 C6 Age 4oC, 3 months 4oC, 3 months 22oC, 3 months 22oC, 3 months 15oC, 2 months 15oC, 2 months Head space 25% 0% 0% 25% 25% 0%  c. Sample group D for tests on age and head space effects (aged pellets prepared for tests at 40oC)  D1 D2 D5 D6 D7 D8 Age 22oC , 3 months 22oC, 3 months 4oC, 3 months 4oC, 3 months 22oC, 2 months 15oC, 2 months Head space 25% 0% 0% 25% 0% 0%  d. Sample group E for age and head space effects (aged pellets prepared for test carried out at 50oC) Sample  E1 E2 E3 Age 22oC , 3 months 4oC, 3 months 15oC, 3 months Head space 25% 0% 0%    5.2. Off-gassing and emission factor Experimental setup and equipment for off-gassing emission measurements have been described in Chapter 2. A 25 ml of gas sample was taken from the test reactor (15 ml for the 2 L container) daily and the gas composition was analyzed by a gas chromatograph, GC-14A, Shimadzu Corporation, Japan). The GC is calibrated for CO, CO2, CH4, H2 O2 and N2 in volume percentage. The gas concentration can be easily influenced by many parameters, such as the amount of pellets, volume of the storage, and void factor in the storage. Gas emission factor (f) was introduced by 81 Kuang et al [33-35] for comparison of their off-gassing experimental results. The emission factor of gas i (fi) which is defined as grams of gas emission per kilogram of pellets,  can be calculated  from gas concentration by the following equation: 0( )i g wt N i P Nt P CV M C f RTM C   . (5.1) where P is the gas pressure, in pa; Ci is the volume fraction (concentration) of gas species i; Vg is the gas volume, which can be calculated by container volume minus the pellets volume (Vg=V-Vpellet), in m3 ; The pellets occupied volume (Vpellet) equals the weight of loaded pellets (Mp) divided by the single pellets density; Mwt is the molecular weight of gas i;  R is the gas constant, 8.314 J/(mol K);  T is the test temperature, in K; MP is the weight of the loaded pellets, in kg. During the test, the total molar number of the gas is reducing everyday due to sampling and reactions. Assuming that the amount of N2 remains constant during off-gassing process as N2 is not generated or assumed, CN0/CNt is introduced as a correction factor to balance the gas volume change. 5.3. Results and discussion 5.3.1. Effects of reactor size and head space on the gas emissions Experiments were carried out  using fresh pellets at 22oC  in three different sizes of reactors, 2 liter reactor (C13, C14), 10 liter reactor (C1, C2) and 45 liter reactor (A1, B1), as shown  in Figure 2.2. Information of pellet samples are listed in Table 5.1a. The head spaces in tests B1, C1 and C14 are 25% and in tests A1, C2 and C13 are 0%. The gases were sampled and analyzed for up to 65 days. The gases detected from these tests are mainly CO2, CO, CH4 and H2. The concentrations and emission factors of those components are plotted in Figure 5.1, Figure 5.2, Figure 5.4 and Figure 5.5. 82 0 10 20 30 40 50 60 70 0.00 0.01 0.02 0.03 C O 2  E m is s io n  f a c to rs  ( g /k g ) Time (days) b 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0     A1, 45L, 0% HS  B1, 45L, 25% HS  C13, 2L, 0% HS  C14, 2L, 25% HS C O 2  C o n c e n tr a ti o n  ( % ) a C1, 10L, 25% HS  C2, 10L, 0% HS  Figure 5.1 Effects of reactor size and head space on (a) CO2 concentration and (b) emission factors for fresh pellets at 22oC. 83 0 10 20 30 40 50 60 70 0.000 0.005 0.010 0.015 0.020 Time (days) C O  E m is s io n  f a c to rs  ( g /k g ) b 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5   A1, 45L, 0% HS  B1, 45L, 25% HS  C13, 2L, 0% HS  C14, 2L, 25% HS C O  C o n c e n tr a ti o n  ( % )  C1, 10L, 25% HS  C2, 10L, 0% HS a  Figure 5.2 The effects of reactor size and head space on CO concentration and emission factors for fresh pellets at 22oC  84 0 10 20 30 40 50 60 -2 0 2 4 6 8 10 12 14 16 18 Time (days)  A1,45L,0% HS  B1,45L,25% HS O x y g e n  c o n c e n tr a ti o n  ( % )  C1,10L,25% HS  C2,10L, 0% HS  C13,2L,0% HS  C14,2L,25% HS  Figure 5.3 The effects of reactor size and head space on O2 depletion for fresh pellets at 22oC. Figure 5.1 shows that the CO2 concentration increased with the reaction time. A maximum CO2 concentration of 2.3% was detected in 45 L reactor with 0% head space. The concentration increased rapidly at the first test run, and reached a steady peak concentration after about 15 days, Figure 5.3 shows the O2 concentration change during the test. O2 concentration drops very fast at the beginning 10 days and reaches to almost zero in 15 to 20 days depending on the headspace and reactor sizes. Because CO2 and CO are the reaction products of the oxidation of the hydrocarbons of the pellets, the oxidation reaction rate is high at the beginning of the run when the oxygen concentration is high, and slows down with  the reaction time when the oxygen is lowered. The reaction stops when all oxygen is consumed. [34, 128]. The production of CO2 can be expressed as a first order chemical reaction [4,6] and its concentration can be expressed by the following first order reaction kinetics [33-35]  ,( ) [1 exp( )]i i iC t C k t   .  (5.2) 85 where Ci(t) is the gas (i) concentration at time t, %; Ci,∞ is the equilibrium concentration, %; ki is the reaction rate constant of the reaction; t is the reaction time, days. The CO2 emission factor (fi) is calculated following Eq. (5.1) and also illustrated in Figure 5.1b. Figure 5.1b shows that similar increasing trends of emission factors with time as for the CO2 concentrations in Figure 5.1a Thus, based on the same assumption, the emission factor (fi) can also be expressed by following first order reaction kinetics  , ( ) [1 exp( )]i i if t f k t   (5.3) where fi(t) is the gas (i) emission factor at time t, g(gas)/kg(pelllets); fi,∞ is the equilibrium emission factor, g(gas)/kg(pellets); ki is the reaction rate constant of the reaction; t is the reaction time, days. All these parameters were calculated and listed in Table 5.2 and Table 5.3. Table 5.2 Summary of emission parameters for CO2 and CO with different reactor sizes and head spaces  CO2  CO   f∞ kf C∞ kC f∞ kf C∞ kC C13 2L, 0% 0.0184 0.1460 1.1351 0.1313 0.0149 0.0605 1.4758 0.0558 C2 10L, 0% 0.0134 0.1185 1.4602 0.1110 0.0081 0.1330 1.3932 0.1227 A1 45L, 0% 0.0213 0.1514 2.3137 0.1478 0.0110 0.2042 1.8801 0.1947 C14 2L, 25% 0.0132 0.1403 1.6931 0.1258 0.0074 0.0941 1.5028 0.0853 C1 10L, 25% 0.0237 0.0665 1.5229 0.0651 0.0137 0.0664 1.3527 0.0865 B1 45L, 25% 0.0273 0.1217 1.8725 0.1169 0.0169 0.1701 1.6283 0.1610  Table 5.3 Summary of emission parameters for CH4 and H2 with different reactor sizes  CH4  H2   f∞ kf C∞ kC f∞ kf C∞ kC C13 2L, 0% 0.0005 0.0340 0.0893 0.0462 0.0003 0.0767 0.9480 0.0540 C2 10L,0% 0.0005 0.0354 0.1503 0.0217 0.0007 0.0530 1.6240 0.0450 A1 45L, 0% 0.0007 0.0838 0.1926 0.0848 0.0006 0.0515 1.4539 0.0504 C14 2L, 25% 0.0003 0.0481 0.0817 0.0308 0.0007 0.0597 0.0816 0.0679 C1 10L,25% 0.0005 0.0367 0.0962 0.0355 0.0007 0.0680 0.9315 0.0665 B1 45L, 25% 0.0008 0.0779 0.1431 0.0736 0.0007 0.0541 1.0952 0.0530  86  Figure 5.1 also shows that the reactor head space has a major effect on gas concentration and emission factor as well. Reactors with 0% head space (A1, C2, C13) have higher CO2 concentrations because of the larger amount of pellets per volume of gas space comparing to the reactor with 25% head space. On the opposite, the emission factors in the reactor with 0% head space (A1, C2, C13) are lower than those reactors with 25% head space because of its smaller gas volume  (i.e. less amount of O2 ). The gas emission factor is more appropriate for reporting and analyzing the off- gassing results because it is based on per unit weight (kg) of wood pellet, eliminating the effect of pellet quantity.  The emission parameters were calculated by fitting data to Eq. (5.2) and Eq. (5.3) using the least square method, with the results being listed in Table 5.2. . Overall, the CO2 emission factor increased with increasing the reactor size for a given head space. The emission factor was the highest in the 45 L reactor (0.027 g/kg in 25% headspace test) and lowest in the 2L reactor, although the difference between 2L and 10L reactors was relatively small. In Figure 5.3, the oxygen concentration for 2L reactor (C13, C14) reached at 1.5% at the equilibrium while other reactors almost reached zero at the equilibrium. This error is caused by the everyday sample procedure, a very small amount of air leaks into the reactor while injection and withdraw of the needle. The volume of the 2L reactor is too small, so that the small amount of air leaked into the reactor caused a system error in 2L reactor. 87 0 10 20 30 40 50 60 70 0.0000 0.0002 0.0004 0.0006 0.0008 Time (days) C H 4  E m is s io n  f a c to rs  ( g /k g ) b 0.00 0.05 0.10 0.15 0.20 0.25 0.30   A1, 45L, 0% HS  B1, 45L, 25% HS a  C13, 2L, 0% HS  C14, 2L, 25% HS C H 4  C o n c e n tr a ti o n  ( % )  C1, 10L, 25% HS  C2, 10L, 0% HS  Figure 5.4 The effect of reactor size and head space on CH4 concentration and emission factors for fresh pellets at 22oC   88 Figure 5.2 shows the effect of reactor size and head space on CO concentration and emission factor for fresh pellets at 22oC. The CO emission profiles were quite similar to those of CO2 emissions, because both were produced from oxidation reactions. The 45 L reactor showed the highest CO concentration, a final concentration of 1.88% was detected in 0% head space test. The CO concentration in 2L reactor and 10L reactor was lower than in 45 L reactor with the same % head space. Emission factors were also calculated and shown in Figure 5.2b. The CO emission factor increased with the reactor size.  The highest emission factor of 0.0169 g/kg(pellets) was measured in the 45L reactor with 25% head space.  The emission factor for the 2L reactor was the lowest, although the difference between 2L and 10L were relatively small. Figure 5.4 shows the effect of reactor size and head space on CH4 concentration and emission factor of fresh pellets at 22oC. A maximum equilibrium CH4 concentration of 0.193% was measured in the 45 liter reactor with 0% head space, which was 10 times smaller than the CO2 and CO concentrations in the same reactor. The concentration of CH4 also increased with the reaction time, rapidly at the beginning and reached a plateau concentration at the end. However, the time to reach the plateau for CH4 was more than 45 days, which was three times longer than CO2 and CO emissions. This is because CH4 is not a product of oxidation reaction which produces CO2 and CO. Therefore, the emission factor of CH4  was less affected by the reactor head space  as shown in Figure 5.4b.       89  0 10 20 30 40 50 60 70 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008  H 2  E m is s io n  f a c to rs  ( g /k g ) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8   A1, 45L, 0% HS  B1, 45L, 25% HS Time (Days)  C13, 2L, 0% HS  C14, 2L, 25% HS H 2  C o n c e n tr a ti o n  ( % )  C1, 10L, 25% HS  C2, 10L, 0% HS  Figure 5.5 The effect of reactor size and head space on H2 concentration and emission factor for fresh pellets at 22oC  90 The reactor size has an important impact on CH4 emission factor as shown in Figure 5.4b. The CH4 emission factor in 45 L reactor was significantly larger than in 10L and 2L reactors.  A maximum plateau emission factor of 0.0008 g/kg(pellets) was reached  in 45L reactor with 25% head space. Larger head space provides more air (oxygen) which promotes the oxygen dependent reactions, such as  the reactions for the formation of CO2 and CO. CH4 emission factor was  not affected by the head space at the beginning of the reaction ; but the emission factor was higher for the reactor with a larger head space as reaction progressed. Besides CO2, CO and CH4, a significant amount of H2 was also detected during the test. Figure 5.5 shows the effect of reactor sizes and head spaces on H2 concentration and emission factor of fresh pellets at 22oC. A maximum H2 concentration of 1.5% was detected in the 10 liter reactor with 0% head space. The concentration of H2 increased steadily   with the reaction time and reached its maximum after more than 60 days, which was even longer than CH4. H2 was produced from non- oxygen dependent reactions; therefore it was less affected by the air to pellets volume ratio (head space) as shown in Figure 5.5b. The results show that the reactor size as well as the head space   did not have a significant effect on H2 concentration and emission factor. A maximum H2 emission factor of 0.0007g/kg(pellets) was found in the 45 liter reactor with a 25% head space. The emission kinetics parameters for CH4 and H2 were also determined by the least square curve- fitting method following Eq. (5.2) and Eq. (5.3),  and listed in Table 5.3. Comparing Table 5.3 to Table 5.2, the reaction constant k of CH4 is much smaller than the reaction constant for CO2 and CO. This also supports that the reactions produced CH4 and H2 (non oxygen dependent reactions) were slower than the oxygen dependent reactions, which produced CO2 and CO. 91 During the experiment, a pressure drop in the sealed reactor was observed. The plastic reactor (10L) was a little distorted and large amount of air was sucked into the reactor when the valve was open after the experiments. This was caused by two reasons. One was the everyday sampling. Withdrawing 25ml gas everyday reduced the gas inside of the reactor. The second was the off- gassing process, For example, in the test with the 45L reactor of 0% head space, about 20% oxygen was consumed when the reaction stopped and yielded concentrations of CO2, CO CH4 and H2 2.3%, 1.88%, 0.193% and 1.5% respectively. During the reaction, the total moles of gases were reduced by 14.4%. 5.3.2. Effect of pellets age Experiments were carried out to study the age effect on the off-gassing emission in this section. The tests were carried out in the 15L and 10L reactors (Figure 2.2) at 22oC. The information of pellet samples are listed Table 5.1b. Pellet age has been defined in Chapter 4 as the days conditioned at a specific temperature. It was also proven that higher conditioning temperature increased the rate pellet aging in Chapter 4. The pellet aging rate   in a cold room at 4oC is really slow, thus it can be assumed to be zero. Therefore, the age of samples C1 and C2, which were conditioned in the cold room for 3 months, is considered as zero; the age of samples C3 and C4 which were conditioned at 22oC for 3 months, is 3 months (22oC), which are the oldest aged samples in the tests; the age of samples C5 and C6 which were conditioned at 15oC for 2 months, is 2 months (15oC) much less than C3 and C4 samples. Figure 5.6, Figure 5.7, Figure 5.9 and Figure 5.10 show the emission factors of CO2, CO, CH4 and H2 for the samples C1-C6 as tested at 22 oC. The emission kinetic parameters were fitted using the least square curve-fitting method following Eq. (5.3) and listed in Table 5.4. 92 0 10 20 30 40 50 60 70 -0.005 0.000 0.005 0.010 0.015 0.020 0.025 0.030   C O 2  E m is s io n  f a c to r (g /k g (p e lle ts ))  C1  C2  Time (days)  C4  C3  C6  C5 25% HS 0% HS  Figure 5.6 CO2 emission factor at 22oC vs. testing days in  10L reactor of 0% and 25% head spaces for pellet samples of different ages.  Figure 5.7 CO emission factors at 22oC vs. testing days in 10L reactor of 0% and 25% head spaces for pellet samples of different ages. 93 0 10 20 30 40 50 60 70 -2 0 2 4 6 8 10 12 14 16 18 20 22  C1  C2 O x y g e n  c o n c e n tr a ti o n  ( % ) Time (days)  C4  C3  C6  25% HS  C5  0% HS  Figure 5.8 O2 depletion at 22oC vs. testing days in 10L reactor of 0% and 25% head spaces for pellet samples of different ages. The age effect on CO2 and CO emission was similar, since both gases were produced from oxidation reactions. Figure 5.6 and Figure 5.7 show that pellets's age has a significant impact on the CO2 and CO emissions. The fresh pellets (C1 and C2) had much higher an emission factor; the 2 month old (15oC) pellets (C5 and C6) had a lower emission factor and the oldest pellets, 3 month old (22oC) (C3 and C4), had the lowest emission factor. The emission factors of CO2 and CO of the fresh pellets (C1 and C2) were significantly affected by the head space, with a much higher emission factor for the reactor with a larger head space, while CO2 and CO emission factors of older pellets were not much affected by the head space. Figure 5.8 shows the O2 concentration decreases during the age effect test. The figure shows that the age of pellets has a very important impact on oxygen depletion. The O2 depleted very fast in the fresh pellets test (C1 and C2) and slower in the 3 month old (15oC) pellets test (C5 and C6) and slowest in the 3 month old (22oC) pellets test (C3 and C4). This phenomenon indicates that the 94 older pellets have lower reaction activities, which we also observed during the self-heating study in Chapter 4. Table 5.4 Summary of emission kinetic parameters for CO2, CO, CH4 and H2 of pellets with different ages  CO2 CO  CH4  H2   f∞ kf f∞ kf f∞ kf f∞ kf C1 fresh, 25% 0.0237 0.0665 0.0137 0.0664 0.0005 0.0367 0.0007 0.6802 C2 fresh,0% 0.0134 0.1185 0.0081 0.1330 0.0005 0.0354 0.0007 0.0530 C3 22oC,0% 0.0118 0.0415 0.0080 0.0192 0.0001 0.0243 0.0003 0.0332 C4 22oC,25% 0.0122 0.0505 0.0130 0.0112 0.0001 0.0129 0.0004 0.0223 C5 15oC 25% 0.0149 0.0854 0.0085 0.0377 0.0003 0.0235 0.0006 0.0305 C6 15oC 0% 0.0127 0.0740 0.0111 0.0260 0.0005 0.0129 0.0007 0.0254  The maximum plateau emission factors, f∞, of CO2 and CO are 0.0237 and 0.0137 g/kg(pellets), respectively, for  C1 sample, the fresh pellets tested with a 25% head space. The plateau emission factor decreased with increasing pellets age at the same head space. Similar trend was observed for the reaction rate constant. Thus, older pellets decompose slower than the fresh pellets and produce less CO2 and CO emissions. This phenomenon has also been observed in Chapter 4 where the oxygen dependent reaction heat was found to decrease when the pellets became older.  95 0 10 20 30 40 50 60 70 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007  C1  C2  C H 4  E m is s io n  f a c to r (g /k g (p e lle ts ))  C4  C3 Time (days) 25% HS 0% HS  C6  C5  Figure 5.9 CH4 emission factor at 22oC vs. testing days in  10L reactor of 0% and 25% head spaces for pellet samples of different ages. 0 10 20 30 40 50 60 70 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008  C1  C2 H 2  E m is s io n  f a c to r (g /k g (p e lle ts ) )  Time (days)  C4  C3 25% HS 0% HS  C6  C5  Figure 5.10 H2 emission factor at 22oC vs. testing days in  10L reactor of 0% and 25% head spaces for pellet samples of different ages. 96  CH4 and H2 emission factors of different aged pellets are listed in Figure 5.9 and Figure 5.10. The CH4 and H2 emissions were also highly affected by the pellet age. The fresh pellet (C1 and C2) had the highest CH4 and H2 emission factors, which are approximately 4 times larger than emission factors of the 3 month old (22oC) pellets (C3 and C4). The CH4 and H2 emission factors show almost identical results for   25% and 0% head space tests in Figure 5.9 and Figure 5.10. The increase in the reactor head space increases the amount of oxygen in the reactor, which could significantly increase the oxygen dependent reactions. Such as those leading to the formation of CO2 and CO; but had less effect on oxygen independent reactions, such as those leading to the formation of CH4 and H2 . 5.3.3. Effect of test temperature The off-gassing emission is very sensitive to the test temperature. Kuang et al. [33] found that the emission of CO2 and CO was  highly dependent on the test temperature in the range of 10 oC to 45oC for fresh pellets. In this section, we tested the off-gassing behavior of two different aged pellets, fresh (C2, D5, E2) and 3 month old (22oC) (C3, D7, E3) pellets at testing temperatures of 22, 40 and 50oC. C2 and C3 samples were tested in 10L reactors while D5, D6, E2 and E3 samples were tested in 45 L reactors.  The emission factors of CO2, CO, CH4 and H2 are shown in Figure 5.11 to Figure 5.15. The emission kinetic parameters were calculated and listed in Table 5.5. 97 0 10 20 30 40 50 60 70 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 50 o C test40 o C test  C2  C3  C O 2  E m is s io n  f a c to r (g /k g (p e lle ts ) ) 22 o C test Fresh pellets Aged pellets Time (days)  D5  D7  E2  E3  Figure 5.11 Effect of  temperature  on CO2 emission factor for fresh and aged pellets 0 10 20 30 40 50 60 70 0.000 0.005 0.010 0.015 0.020  C2  C3  C O  E m is s io n  f a c to r (g /k g (p e lle ts ) )  D5  D7 50 o C test40 o C test  E2  E3 Time (days) Fresh pellets Aged pellets 22 o C test  Figure 5.12 Effect of temperature on CO emission factor for fresh and aged pellets 98 Figure 5.11 and Figure 5.12 show the emission factors of CO2 and CO at testing temperatures of 22, 40 and 50oC. The ages of pellets showed less effect on the emission of CO2 and CO at 40 oC and 50oC, comparing to at 22oC, especially for the CO emission. For the test at 22oC, the aged pellets C3 had a much less CO emission than the fresh pellets C2.   However the aged pellets tested at 40oC and 50oC (D7 and E3) had almost the same emissions as the fresh pellets (D5 and E2). Same trends can be found in emissions of CH4 and H2 in Figure 5.14 and Figure 5.15. Figure 5.13 shows that the O2 depletion speeds in fresh pellets test were faster than aged pellets test. At higher temperature (40oC and 50oC), the O2 depletion speeds difference between fresh pellets and aged pellets were much smaller than at 22oC test. In another word, the age of pellets is an important factor that could impact the gas emission substantially; however, it has less impact when the test temperature is higher. This observation is consistent with the effect of pellet ages on the self-heating of pellets during storage (see Chapter 4). 0 10 20 30 40 50 60 70 0 5 10 15 20 25 50 o C test40 o C test O x y g e n  c o n c e n tr a ti o n  ( % ) Time (days)  C2  C3 22 o C test  D5  D7  E2  E3 Fresh pellets Aged pellets  Figure 5.13 Effect of temperature on O2 depletion for fresh and aged pellets   99 Table 5.5 Summary of emission kinetic parameters for CO2, CO, CH4 and H2 of pellets at different tested temperature     CO2 CO   CH4   H2   f∞ kf C∞ kC f∞ kf C∞ kC E2  0.0565 0.6230 0.0130 1.0723 0.0049 0.0329 0.0008 0.5802 E3  0.0589 0.4026 0.0136 0.5310 0.0030 0.0840 0.0008 0.3692 D5  0.0570 0.1492 0.0090 0.2596 0.0013 0.8370 0.0008 0.1296 D7  0.5916 0.2006 0.0084 0.2064 0.0012 0.0799 0.0006 0.1481 C2  0.0134 0.1185 0.0081 0.1330 0.0005 0.0354 0.0007 0.0530 C3  0.0118 0.0415 0.0069 0.0221 0.0001 0.0243 0.0003 0.0332  Figure 5.14 and Figure 5.15 show the emission factors of CH4 and H2 at testing temperatures of 22, 40 and 50oC using two different aged pellets. The emission factors of CH4 and H2 increased with increasing the test temperature. The effect of pellet age on the emission factor was reduced at higher test temperatures. Even though the high temperature increases the plateau emission factor and the reaction rate constant substantially, the reaction rate constants for CH4 and H2 were still much smaller than those for CO2 and CO.  100 0 10 20 30 40 50 60 70 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 0.0016  C2  C3  C H 4  E m is s io n  f a c to r  ( g /k g (p e lle ts ) ) Time (days) 50 o C test40 o C test  D5  D7 Fresh pellets Aged pellets 22 o C test  E2  E3  Figure 5.14 Effect of temperature on CH4 emission factor for fresh and aged pellets. Temperature is added in the legend in the following figures according to dr. Bi's comment 0 10 20 30 40 50 60 70 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0.0014  C2  C3 H 2  E m is s io n  f a c to r (g /k g (p e lle ts ) ) Time (days)  D5  D7 50 o C test40 o C test  E2  E3 22 o C test Fresh pellets Aged pellets  Figure 5.15 Effect of temperature on H2 emission factor for fresh and aged pellets 101 5.4. Conclusions Experiments were carried out to study the gas emission during bulk storage of wood pellets. The reactor size and head space, pellets age and test temperature on off-gassing were investigated in this chapter and the following conclusions are obtained 1. The off-gassing from bulk storage of wood pellets mainly consists of CO2, CO, CH4 and H2. The gas emissions approximately follow a first order reaction kinetics with maximum concentrations of 2.3%, 1.88%, 0.193% and 1.5%, respectively, and plateau emission factors of 0.0273, 0.0167, 0.0008 and 0.0007 g/kg(pellets), respectively. 2. The emission factors of CO2, CO and CH4 increase when the reactor size is larger and the emission factor of H2 is less affected by the reactor size. 3. Pellets age has significant impacts on CO2, CO, CH4 and H2 emissions. Emission factors of all gases decrease with increasing the age of pellets. However, the effect of pellets age decreases as the test temperature increases. 102  6.  Determination of self-heating using lab scale experiments1 This chapter presents the results of self-heating kinetic parameters at high temperatures determined by the lab scale experiments. 6.1. Materials Softwood pellets produced in British Columbia (BC) were used in this research. Freshly made pellets with designation A, B, C and E in this research were obtained from 2 different wood pellet plants in BC. They were made of whitewood mostly pine.  Pellet A and pellet C were manufactured by Princeton Co-Generation Corp. (Princeton, BC, Canada). Pellets B and E were manufactured by Premium Pellet Ltd. (Vanderhoof, BC, Canada). Pellets were on average 6.2 mm in diameter and 10 to 25 mm in length. The pellets were made from sawdust and shavings, by-products of milling spruce, pine and fir for lumber. Upon their arrival at the Lab, the pellets were dried for 24 hours at 103oC to remove the moisture before each test. Experiments were carried out to investigate self-heating in 4 different sizes of reactors placed in a laboratory oven (Thermal scientific, 1.7 cu. ft). The experimental setup is described in section 2.4. The containers were filled up with wood pellets (weight was recorded) and placed  inside the oven maintained at a constant temperature between 120 and 200oC, with temperatures being recorded by a data acquisition system. According to FK method, experiments were carried out to obtain the critical oven temperature by a number of trial and error tests. First,  the reactor (filled with fresh wood pellets) was placed into the oven, then the oven was set to a target temperature and maintained at this temperature for at least 24 hours. If there was a thermal runaway, the test was repeated using a lower oven temperature. Otherwise, the test was continued using a higher oven temperature, until a minimum oven temperature is obtained and this oven temperature is defined as the critical oven temperature triggering the pellet temperature runaway for this particular reactor.  1 A version of this chapter will be submitted for publication. Guo, W., Lim, CJ., Bi, X., Sokhansanj, S., Melin, ,S.. Determination of self-heating generation rate of wood pellets in storage using a lab scale experiment.. 103 The experimental results were analyzed using FK method, CP method and numerical curve fitting method in order to evaluate the kinetic parameters (E and ΔrhA). A two-dimensional self-heating model was developed and solved numerically using using Matlab toolbox. 6.2. Thermal runaway and spontaneous combustion The thermal runaway and spontaneous combustion phenomena were investigated in a top open cylindrical aluminum reactor of 417 mm in height and 356 mm in diameter. The reactor was filled with pellets and maintained in an oven with a constant temperature. Four K-type thermocouples of 3 mm in diameter were inserted in the bulk pellets for monitoring the temperature during the experiment. The thermocouples (T1, T2, T3, T4) are inserted from the top with their tips being maintained at the same axial level 200 mm above the bottom and their distances to the center were 116mm, 58mm, 0mm and 88 mm, respectively. The experimental diagram is illustrated in Figure 2-5. Figure 6.1 shows a set of measured temperature profiles. When the oven temperature was set at 90oC, T1 increased faster than others, and reached a plateau temperature of 84oC after 24 hours, as shown in Figure 6.1a. In other words, the temperature difference between the reactor wall and T1(r=116mm) was very large while the difference between T1 (r=116mm) and the center was small. No self-heating was detected in this test. When the oven temperature was set at 105oC as shown in Figure 6.1b, T1 showed the fastest increase with time because T1 was the closest to the reactor wall. T4 (r=88mm) was lower than T1 because its location was further away from the heating source than T1. T3 was located in the center, which was expected to have the lowest temperature if no self- heating. However, T2 and T3 were much higher than T4. This means that heat generated from wood pellets itself started to become detectable at this temperature and to cause the center temperature increase at this test segment. When the oven temperature was set at 125oC as shown in Figure 6.1c, the heat generation from pellets became more obvious and caused the center temperature (T3) to cross  T2 (r=58 mm) at around 115oC. The temperature at this point is defined as the crossing point (CP) temperature which is an important point to be used to identify the development  of self-heating. When the oven temperature increased to 132oC in Figure 6.1d,  self- heating was even more significant causing the center temperature (T3) to cross   T1 (r=116 mm) and to become the highest temperature in the reactor. It is noted that the highest temperature in the 104 reactor is still lower than the oven temperature, which means that self-heating at this test segment is not sufficient to cause thermal runway. 0 200 400 600 800 1000 1200 1400 1600 20 30 40 50 60 70 80 90 100 Time, (min) T em pe ra tu re , ( o C ) SI11-18-2.xls   Toven T1 T2 T3 T4  a 0 500 1000 1500 80 85 90 95 100 105 110 Time, (min) T em pe ra tu re , ( o C ) SI11-19-2.xls   Toven T1 T2 T3 T4  b 0 500 1000 1500 2000 2500 3000 95 100 105 110 115 120 125 Time, (min) T em pe ra tu re , ( o C ) SI11-20-2.xls   Toven T1 T2 T3 T4  c 0 500 1000 1500 2000 2500 3000 116 118 120 122 124 126 128 130 132 134 Time, (min) T em pe ra tu re , ( o C ) SI11-21-2.xls   Toven T1 T2 T3 T4  d 0 500 1000 1500 2000 2500 3000 3500 120 140 160 180 200 220 240 Time, (min) T em pe ra tu re , ( o C ) SI11-24-2.xls   Toven T1 T2 T3 T4  e 0 1000 2000 3000 4000 5000 6000 0 50 100 150 200 250 Time, (min) T em pe ra tu re , ( o C ) SI11-26-2.xls   Toven T1 T2 T3 T4  f Figure 6.1 A set of measured temperature profiles from the self-heating reactor. 105 the thermocouples (T1, T2, T3, T4) are inserted from the top with their tips being maintained at the same axial level 200 mm above the bottom and their distances to the center being 116 mm, 58 mm, 0 mm and 88 mm, respectively. The experiment setup is illustrated in Figure 2-5. Notes: the pellets are pre-dried When the oven temperature increased to 143oC in Figure 6.1e, temperatures of all locations ran away from the oven temperature and increased dramatically. This oven temperature was identified as the critical temperature at which the thermal runaway happened. At 2500 minutes, all temperatures went up from temperature 160oC to230oC very quickly with the center temperature (T3) being the highest, T2 (r=58 mm) being the second highest and followed by T3 (r=88 mm) and T4 (r=116 mm). This phenomenon is identified as the self-combustion (ignition). After the ignition, the oven turned off and, as the temperature profiles in Figure 6.1f show, the combustion continued and temperatures fluctuated between 220oC and 150oC for 4 days with no tendency of cooling down. In this case, the heating source is only from the wood pellet itself, as evidenced by the temperature gradients from the center toward the outside wall with the center (T3) temperature being the highest and T4 (closest to the wall) being the lowest. After the reactor center had been combusting for 4 days, the reactor was removed from the oven, with the appearance of the pellets shown in Figure 6.2. It is interesting to note that there was no smoke or smell coming out of the reactor during the combustion.  The pellets on the surface were still in the original form and no sign of burning at all except some moisture gathered near the center thermocouple, as shown in Figure 6.2a. However, the pellets about 15 cm below the surface were charring and emitting heat, with no light or smoke during initial stages of exposure to outside air, as shown in Figure 6.2b. This phenomenon has been discussed in Chapter 1 as smoldering combustion (section 1.2.4). As shown in Figure 6.2b, the smoldering started from the center and gradually spread out. With limited oxygen supply, the cellulosic materials underwent   pyrolysis   that  produced active char and large amount of heat, instead of the gas oxidation combustion, commonly known as burning or flaming. However, the smoldering phenomenon can easily develop   to a rapid flame spreading combustion once the air flow increases or other condition changes which might increases the oxidation rate [64, 65]. As we observed, the pellets started smoking after exposed in air for a couple of minutes and can easily develop to flame combustion. . 106  a. The surface of the burning pellets.  b. Pellets are burning 15 cm below the surface Figure 6.2 Pellets surface after self-heating experiment.  107  6.3. Crossing point and critical oven temperature In this research two important temperatures are defined, crossing point and critical oven temperature. These two temperatures are critical in determining whether a bulk pellet will reaches the point of smoldering and spontaneous combustion. From chapter 4, we learned that self-heating rate is small at room temperature and increases at temperatures higher than 40oC. At lower temperatures or in a small container, released heat can dissipate quickly from the bulk pellets so that a stationary state will be reached and no crossing point appears in the temperature profile such as shown in Figure 6.1a. The rate of self-heating increases at higher temperatures because the heat cannot dissipate quickly. This results in the center temperature surpassing the off-center temperature. The crossing point is reached when the center temperature crosses the off center temperature. The crossing point is an important parameter to indicate whether self-heating becomes significant enough to cause temperature increase inside the bulk material. Figures 6.1b-d show the evolution of temperatures inside the bulk especially when the crossing point was reached at 96oC (T2 and T3) in Figure 6.1b, 115oC (T2 and T3) in Figure 6.1c, 132oC (T1 and T3) in Figure 6.1d and 144oC (T1, T2, and T3) in Figure 6.1e. As Figures 6.1 shows, reaching the crossing point depended upon the amount of heat input to the container. When self-heating was high enough, a thermal runaway appeared in the temperature profiles and the temperature shoot up to 200oC, which indicated the occurrence of spontaneous smoldering combustion in the bulk wood pellets. The critical oven temperature is defined as the minimum oven temperature that leads to a thermal runaway in the temperature profile and causes the bulk pellets temperature rise to over 200oC. Following the same procedure, the crossing point and thermal runaway temperatures are determined for different size of reactors. The results are listed in Table 6.1 and some sample results are shown in Figure 6.3. The reactors are top open cylindrical reactors, with two thermocouples in each of them, T1 (Center) and T2 (off-center). The reactor sizes are Large (H=224 mm, D=216 mm), Medium (H=188 mm, D=164 mm), Small (H=114 mm, D=105 mm). The details of the reactor setups are described in section 2.4. The pellets are pre-dried before a self-heating experiment. 108  Table 6.1 The steady state center temperature for 3 sizes containers at different oven temperatures (Pellet E). Oven temperature Large reactor  Medium reactor Small reactor FS CP Max T FS CP Max T FS CP Max T 100oC SS  100 oC SS  100 oC SS  100 oC 130oC SS yes 132 oC SS  130 oC SS  130 oC 140oC SS yes 144 oC SS yes 141 oC SS  140 oC 145 oC    SS yes 148 oC 150oC SS  162 oC SS  153 oC SS  150 oC 160oC TR  196 oC SS  169 oC SS yes 162 oC 170oC TR  - SS  178 oC  SS  179 oC 180 oC TR  - TR  195 oC  SS  188 oC 191 oC  TR  - TR  - TR  200 oC Notes: SS stands for steady state being reached; TR stands for a thermal runaway being observed in the experiment; Max T stands for the highest temperature observed in the experiments.  A crossing point appears in the large reactor when the oven temperature is 130oC; and after the crossing point, a stationary sate is reached with a maximum temperature of 132oC observed in the center of the reactor as shown in Figure 6.3a. A thermal runaway is observed in the large reactor when oven temperature is 160oC as shown in Figure 6.3d.  The temperature reached 196oC within 2 days and kept increasing after the oven was turned off. Therefore, the critical oven temperature for the large reactor is 160oC. The crossing point for the medium reactor was observed when the oven temperature was 140oC as shown in Figure 6.3b. A steady state was reached with a maximum temperature of 141oC being observed after the crossing point at this oven temperature. Both the crossing point temperature and critical oven temperature increased when the reactor size decreased. A crossing point was observed in the small reactor when the oven temperature was at 160oC.  The thermal runaway was observed in the medium reactor and small reactor when the oven temperature was at 180oC and 191oC, respectively. Because of the low effective thermal conductivity of bulk wood pellets, the larger a reactors is, the more difficult for the heat to dissipate from the bulk pellets to outside and, thus, the lower is the temperature of the crossing point and less chance for thermal runaway. 109 While the crossing point can vary depending on the initial temperature of the pellets, the critical oven temperature is unique for a given reactor since it is the minimum oven temperature or environmental temperature  which leads to the thermal runway. Experiments were carried out to determine the unique critical oven temperature of three reactors by using a trial and error method. The resulted critical oven temperatures are 160 oC for the large reactor, 173 oC for the medium reactor and 191 oC for the small reactor. Those temperatures are used in Frank-kamenetskii method to determine the self-heating kinetics.  0 100 200 300 400 500 600 700 800 900 1000 100 105 110 115 120 125 130 135 Time, (min) T e m p e ra tu re , (o C )   Toven L-T1 L-T2 M-T1 M-T2 S-T1 S-T2  0 200 400 600 800 1000 1200 128 130 132 134 136 138 140 142 144 146 Time, (min) T e m p e ra tu re , (o C )   Toven L-T1 L-T2 M-T1 M-T2 S-T1 S-T2  a. oven temperature 130oC b. oven temperature 140oC 0 500 1000 1500 2000 2500 135 140 145 150 155 160 165 Time, (min) T e m p e ra tu re , (o C )   Toven L-T1 L-T2 M-T1 M-T2 S-T1 S-T2  0 200 400 600 800 1000 1200 1400 1600 140 150 160 170 180 190 200 Time, (min) T e m p e ra tu re , (o C )   Toven L-T1 L-T2 M-T1 M-T2 S-T1 S-T2  c.  oven temperature 150 oC d.  oven temperature 160 oC Figure 6.3 Typical temperature profiles of three different sizes of reactors with increasing the oven temperature.  110 6.4. Determination of self-heating kinetic parameters The self-heating kinetic parameters are determined using three methods in this study. As introduced in Chapter 1, Frank-kamenetskii (FK) method and crossing point (CP) method are two widely used methods to determine the self-heating kinetic parameters. During the self-heating before the thermal runaway, the energy balance within the reactor can be expressed as  2 2 exp( / )P r T T C h A E RT t x              (6.1) where Cp is the specific heat capacity of the pellets, J/(kg.K),  λ is the thermal conductivity of pellets , W/(m.K);  t is time, s; ∆rh is the reaction enthalpy of  self-heating  , J/kg; ρ is bulk density of the pellets, kg/m3; A is the dimensionless pre-exponential factor; E is the activation energy of self-heating reaction, J/mol; R is the gas constant, =8.314 J/(mol K) ;  T is the temperature of the pellets, K. Equation 6.1 assumes  that heat flow is in a single direction x from center towards outside . 6.4.1. Frank-kamenetskii method Frank-kamenetskii (FK) method [102] was developed based on one dimensional slab. Under the steady state, the transient temperature term is eliminated (∂T/∂t=0) and the energy balance equation becomes  2 2 exp( / )r T h A E RT x           (6.2) By introducing two dimensionless factors, 20 0( / )( )E RT T T   , / cr r  , Eq. (6.2) is converted to a dimensionless form,  2 2 e         (6.3) in which the dimensionless parameter δ is defined as 111  0 2 2 0 E RTrhAE r e RT       (6.4) For values of δ where equation Eq. (6.4) has a solution, a stationary solution is possible and the stationary temperature distribution can be obtained. For those values of δ where Eq. 6.3 has no solution, the stationary state cannot be reached and an explosion (thermal runaway in this situation) will take place. Therefore, the maximum δ value which gives a stationary solution (δC) describes the critical condition for the thermal runaway. The critical δc value has been found for different geometry, which can be found in Table 1.5 [101, 102]. The height and diameter ratio (H/D) of the three cylindrical reactors are 1.03, 1.14, 1.08 for larger, medium, and small reactor, respectively. Therefore, it is reasonable to use the critical condition δC for a short cylinder of diameter equal to height, which is found to be 2.844 from Table 1.5. T0, the critical oven temperature (the minimum oven temperature leads to the thermal runaway) has been experimentally determined for each reactor in section 6.3. The effective thermal conductivity λ of regular white wood pellets has been determined in   chapter 3 as  (0.219 0.01 ) (1 ) 0.027eff M         (6.5) where λeff is the effective thermal conductivity, in W/(m K); M is the moisture content, % wet based which is assumed to be 0 in this study since the pellets are pre-dried By applying a natural logarithm, Eq. (6.4) become  0 0ln 2ln( / ) ln( ) / rE hAT r E RT R        .  (6.6) The values of E and ∆rhA are obtained from a linear regression by plotting 0ln( ) 2ln( / )C T r   as y coordinate and 01/ T as x coordinates. Figure 6.4 shows the linear regression of pellets A and Pellets E using FK method. 112 0.00217 0.00224 0.00231 0.00238 17.5 18.0 18.5 19.0 19.5  Pellet E ln ( C )+ 2 ln (T /r ) 1/T 0 y= -9849x+40.40 R 2 =0.993  0.00217 0.00224 0.00231 0.00238 17.5 18.0 18.5 19.0 19.5  Pellet A ln ( C )+ 2 ln (T /r ) 1/T 0 y= -9647x+40.43 R 2 =0.995  a. Pellet E b. Pellet A Figure 6.4 Linear regression of 0ln( ) 2ln( / )C T r   vs 01/ T  for pellets A and Pellets E using FK method. Note: T0, the critical oven temperature, K;  The regression results are listed in Table 6.2. The activation energy E was determined  to be 81.9 and 80.2 kJ/mol for pellets E and pellets A, respectively; the ∆rhA values were determined to be 5.45×106 and 5.71×106 kJ/(kg s) for pellets E and pellets A respectively. The regression results show a high coefficient of determination with R2 = 0.99, which indicates the reliability of this method. 6.4.2. Crossing point method Crossing point (CP) method was developed based on the transient solution of Eq. (6.1). According to the crossing point method, when two temperatures inside of a container become identical (crossing point), it is assumed by definition that there will be no conductive heat transfer between these two points, and equation (6.1) becomes  exp( / )P r P P T C h A E RT t          (6.7) where p indicates the crossing point; TP is the crossing point temperature, K; / PT t  is the temperature increase rate at the crossing point, K/s; The value TP and / PT t   are experimentally determined. By applying a natural logarithm, Eq. (6.7) becomes 113  ln( ) ln( )r P P P hAT E t C RT      (6.8) Following equation (6.8), activation energy E is obtained from a linear regression by plotting ln( / ) P dT dt  vs 1/ PT , The specific heat CP of regular white wood pellets has been determined in our previous study[125] and combine with the literature values [85, 86] at the high temperature, the specific heat at higher temperature can be correlated by P,1C = (1- M /100)(0.00546T -0.524) + 0.042M (6.9)  where CP is specific heat, kJ/(kg.K); M is moisture content in weight percentage which is assumed to be 0 in this study since the pellets are pre-dried; T is  temperature, K; Figure 6.5 shows the linear regression of pellets A and Pellets E using CP method. 0.0021 0.0022 0.0023 0.0024 0.0025 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 ln (d T /d t) /l n (K /s )  Pellet E 1/Tp (1/K) y= -8123x+11.77 R 2 =0.537   a. Pellets E b. Pellets A Figure 6.5 Linear regression of ln( / ) P dT dt  vs 1/ PT   for pellets A and Pellets E using CP method. The regression results are listed in Table 6.2. The activation energy E was estimated  to be 67.5 and 54.6 kJ/mol for pellets E and pellets A, respectively. The ∆rhA values were determined to be 1.81×106 and 7.31×106 kJ/(kg s) for pellets E and pellets A, respectively. These values are smaller than the values from FK method and also have a much lower coefficient of determination (R2), 0.537-0.659. The lower coefficients of determination indicate that the CP method  is less reliable than FK method. This is caused by the different assumptions in CP method and FK method. In CP 114 method, it is assumed that there is no heat conduction between the two points when they reach identical temperatures (the crossing point), which may not true in the real process. This inaccurate assumption may have caused the scatter of the data during the regression in Figure 6.5 and also caused the lower coefficients of determination and lower parameter values in Table 6.2. Table 6.2 E and ∆rhA values determined from different methods.  FK method CP method Curve-fitting Pellet E Pellet A Pellet E Pellet A Pellet E* Pellet A* E, kJ/mol 81.9 80.2 67.5 54.6 78.7 77.4 ∆rhA kJ/(kg s)  5.45*106 5.71*106 1.81*105 7.31*103 2.88*106 2.82*106 R2 0.993 0.995 0.537 0.659 0.985 0.95 *   single step curve fitting, temperatures range from 100 to 180 oC (oven temperature 170 oC ).  6.4.3. Numerical curve fitting method Both the FK method and CP method used a simplified form of energy balance equation Eq. (6.1) and only one specific point in the temperature profile. Frank-Kamenetskii method used the steady state of the energy equation and the critical oven temperature of different sizes reactors, while crossing point method used the transient energy equation (eliminated the heat conduction) at the crossing point. In this study, we also developed a numerical curve-fitting method to determine the kinetic parameters using un simplified energy balance equation Eq. (6.1) and the entire set of temperature profile data, instead of just one point. Eq. (6.1) is solved numerically using Matlab by assuing the reactor is an axi-symmetric and one dimensional cylinder. Curve-fitting was conducted by minimizing the difference between the measured temperature data and the numerical solution. The kinetic parameters (E and ΔrhA) are thus determined by minimizing the residuals between the numerical solution and the experimental temperature profiles. Boundary conditions and assumptions used in the fitting are given below. Assumptions: 1) Bulk pellets can be treated as homogeneous materials with effective thermal and physical properties (Cp, λ, ρ). 2) Initial condition t=0, T=Ti (initial temperature). 115 3) No air movement inside the container. 4) Infinitely large thermal conductivity of the reactor wall (Twall=Toven). Table 6.3 lists the estimated parameters for each reactor in the temperature range of 100 to 170oC. The activation energy E was determined to be 78.3, 78.7 and 80 kJ/mol for large, medium and small reactors, respectively. The ∆rhA values were determined to be 3.39×10 6, 2.69×106 and 9.33×105 kJ/(kg s) for large, medium and small reactor, respectively. The E values from medium and large reactors are close, while the E value from the small reactor is slightly higher. This is probably because the heat dissipates very quickly in the small reactor and the self-heating become more difficult be sustained so that the observed activation energy becomes higher. Because of the similar results, only one reactor (medium size or large size) is needed for using the curve-fitting method to determine the kinetic parameters. As mentioned in the methodology section, the process of determining the critical oven temperature requires large number of trials in order to obtain results for different reactor sizes. Compared to the F-K method, curve-fitting method is less time consuming with a similar accuracy. Table 6.3 E and ∆rhA values determined for different size reactors (pellet E)  Size of container reactor Large * Medium* Small* E, kJ/mol 78.3 78.7 80 ∆rhA, kJ/(kg s)  3.39E+06 2.69E+06 9.33E+05  * Curve fitting in the temperature range of 100-170 oC 116  Table 6.4 shows the kinetic parameters (E and ∆rhA values) of pellets from different manufacturers. The results show a small, but negligible difference between different pellet samples. After tested at different oven temperatures, reactor sizes and pellet samples, a mean activation energy of 78.7±0.8 kJ/mol and a mean ∆rhA of (4.2±2.5)*10 6 kJ/(kgs) are obtained. The activation energy from this study is consistent with the reported values of 76-93 kJ/mol[90] at 100-200 oC and 60 to 80 kJ/mol[26] at 50-80 oC for Swedish wood pellets. Table 6.4 E and ∆rhA values for different types of pellets  Pellet A Pellets B Pellets C Pellets E Mean E, kJ/mol 77.4 77.7 79.6 78.7 78.7±0.8 ∆rhA, kJ/(kg s)  2.82E+06 7.41E+06 9.12E+06 2.88E+06 (4.22±2.5)*106 R2 0.95 0.92 0.95 0.98  The temperature profiles predicted using the mean activation energy of 78.7±0.8 kJ/mol and ∆rhA value of (4.22±2.5)*106 kJ/(kg s) are plotted and compared with the experimental temperature profiles in Figure 6.6. For most of experiments, the average values of E and ∆rhA can fit the data very well, especially in predicting the thermal runaway at a constant oven temperature. However, the predicted thermal runaway sometimes shows up earlier or later slightly compared to the experimental results, such as in Figure 6.6b-d. This is because we assumed a zero order reaction with the reaction rate (A) only depending on the temperature and neglected the depletion of the reaction source and the effect of airflow from natural convection. However, in the experiment, the thermal runaway starts with a slow smoldering reaction and accelerates much more slowly than the prediction.  117  0 100 200 300 400 500 600 700 800 900 1000 0 50 100 150 200 250 300 Time, (min) T e m p e ra tu re , (o C )   Pellet A Pellet A Pellet C Pellet C Numerical Numerical  a. Pellets A and C in medium reactor at oven temperature 170 oC 0 100 200 300 400 500 600 700 80 100 120 140 160 180 200 220 240 260 280 300 Time, (min) T e m p e ra tu re , (o C )   Pellet B Pellet B Numerical Numerical  b. Pellet B in medium reactor at oven temperature 170 oC (initial temperature 100 oC) 0 100 200 300 400 500 600 700 800 900 1000 0 50 100 150 200 250 300 Time, (min) T e m p e ra tu re , (o C )   Pellet E Pellet E Numerical Numerical  c. Pellet E in large reactor at oven temperature 170 oC 0 1000 2000 3000 4000 5000 100 120 140 160 180 200 220 240 260 280 300 Time, (min) T e m p e ra tu re , (o C )   Pellet E Pellet E Numerical Numerical  d. Pellet E in large reactor at oven temperatures of 130, 140, 150 and 160oC Figure 6.6 Temperature profiles predicted using mean activation energy value of 78.7±0.8 kJ/mol and ∆rhA value of (4.22±2.5)*106 kJ/(kg s) 118  6.5. Prediction of the critical ambient temperature Critical ambient temperature is the minimum ambient temperature that will lead to a thermal runaway during the pellet storage. This critical ambient temperature is the same as critical oven temperature. It is important to know the critical ambient temperature for a specific storage container during wood pellets storage in order to prevent the thermal runaway and the spontaneous combustion. Before numerical method was well developed, F-K method was often used as a tool for predicting the critical ambient temperature for a particular reactor. From the linear regression (Figure 6.4) of F-K method, the average critical temperature of wood pellets for different size reactors can be represented as 0 0ln 2ln( / ) 40.4 9748/T r T                                    .(6.10) Following this equation, the relation between the predicted critical ambient temperature (T0) and the container radius is obtained and shown in Figure 6.8. The numerical model also provides a method to predict the relationship between the critical temperature and the reactor size. The numerical model was developed using an axi-symmetric two- dimensional solution of the energy balance equation, with the assumptions and boundary conditions shown in Figure 6.7. In the prediction, the thermal runaway was considered to occur when the center temperature reaches 200 oC and the critical ambient temperature was defined as the environment temperature corresponding to the thermal runaway. The initial temperature of the bulk pellets was set at 10oC below the ambient temperature and the total calculation time was set as n=13×R days (R is the reactor radius, m). The prediction results are plotted in Figure 6.8  119 R L = 2 R            B.C Twall=Ta, at z=L; Twall=Ta, at z=0; Twall=Ta, at r=R; əT/ər=0, at r=0; z r 2 2 1 [ ( ) ] exp( )P r T T T E C r h A t r r r z RT                  Equation:  Figure 6.7 Boundary conditions of the container (filled with wood pellets) Assumptions: heat loss on the container wall is only by conduction of the non-flowing air; the container wall temperature equals the ambient temperature. Note: r is the axial distance to the axis; z is the vertical distance from the bottom. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 0 50 100 150 200 250  Wood fibre  FK method prediction  Experiment data  2D Numerical prediction C ri ti c a l T e m p e ra tu re  ( o C ) Radius (m)  Figure 6.8 Predicted relationships between critical temperature and container radius Note: wood fiber data are experimental data for an octagonal pile of wood fibre [101]; the radius for octagonal pile is defined as the half shortest distance between opposit faces. FK method prediction is obtained from Eq. (6.10); Experimental data are the critial temperature of wood pellets; 120 Figure 6.8 shows the comparison of the critical temperatures predicted by the F-K method and  the 2D numerical method with the experimental data and data reported for wood fibers [101].  The simulation and experimental data show the same trend that the critical temperature that leading to thermal runaway is lowered when the container size (diameter) becomes larger. The critical temperature predicted using the 2D numerical method shows an excellent agreement with the measured critical temperature of wood pellets from this study (experiment data), even though it is lower than the measured critical temperature of wood fiber. The predicted critical temperature from the FK method is much lower than the results from the numerical solution. According to F-K method, 30 oC could be the critical temperature for ignition when the container radius is larger than 5 meters. Since F-K method neglects the transient term and just considers the stationary condition, the numerical solution method is expected to be more reliable for the prediction of critical temperature. The numerical prediction shows that when the radius is larger than 5 meters, the critical temperature is around 80 oC. The model was built based on the experimental results in the temperature range of 100-200 oC and the materials were pre-dried.  The heat release by the biological and chemical reactions happening in the temperature range of 30 oC to 80 oC [26, 129] have not been  taken into consideration in the model. 121  6.6. Conclusions Laboratory experiments were carried out to investigate the self-heating phenomenon at high temperatures (100oC to 200oC). When the ambient temperature is high enough, a thermal runaway appeared in the center of stored wood pellets. When the temperature reached 200oC, the heat generated from the pyrolysis reaction became significant enough to cause an invisible smoldering combustion within the bulk wood pellets. In order to prevent  thermal runaway, this study has evaluated the kinetic parameters (E and ΔrhA) using three methods, F-K method, CP method and numerical curve fitting method for four types of wood pellets. The critical ambient temperature is predicted using both the FK method and an axi-symmetric two-dimensional numerical model.  The following conclusions are obtained 1.The critical ambient temperatures were found to be 135 oC, 160 oC, 173 oC and 191 oC for reactor sizes of radius 0.278, 0.108, 0.082 and 0.053 m, respectively. 2.Both the FK (Frank-Kamenetskii's) method and the curve-fitting method determined the kinetic parameters with a high coefficient of determination. The crossing point method results were too scattered with a much lower coefficient of determination than F-K and numerical methods. Mean values for activation energy (E) of 78.7±0.8 kJ/mol and self- heating rate constant (∆rhA) of (4.22±2.5)*10 6 kJ/(kg s) were obtained for four white pellets samples from different pellet producers in British Columbia. 3.The two-dimensional numerical model shows good agreement with the wood pellets experimental data and more reliable than that predicted by the FK method. The numerical prediction shows that when the radius is larger than 5 m, the critical temperature can be as low as 80 oC.    122 7. Simulation of self-heating in wood pellets silo In this chapter, a self-heating model was developed based on mass and energy transfer equations for a wood pellet silo taking advantage of the thermal properties and the self-heating kinetic parameters collected from the experiments. This model was then used to predict the temperature profiles during the self-heating process, and to identify the thermal runaway condition and time to reach ignition. The effects of wall insulation, outside wind speed, container dimension and inside ventilation were investigated. At the end, the model was verified by comparing the simulation results with temperature profiles from industrial silo, pilot scale experiments and lab scale experiments. 7.1. Self-heating in commercial wood pellets storage silo Self-heating causes temperature increase especially in wood pellets storage silos. Self-heating and accidental fires are observed all over the world. In the previous chapters, we have experimentally measured major parameters affecting the self-heating process, including thermal conductivity, specific heat capacity, and heat generation at both low temperatures and high temperatures. Using those parameters, a two dimensional axi-symmetric model is developed to simulate the self-heating process in wood pellets silos and predict conditions at which a thermal runaway situation might happen. For the purpose of verification, the self-heating data from an industrial silo are collected to compare with the simulation results. Fibreco Inc. is a port bulk handling facility in North Vancouver, British Columbia. The Company has been handling wood pellets for shipment to European markets since 2005. Wood pellets are received by train from BC interior plants, unloaded, transported and stored in six wood pellet storage silos at the site. Each silo has a storage capacity of 4500 metric tons in dimension of 21.9 meters in diameter and 23.2 meters in total height. The pellets are loaded into the silo from the top using tripper-style conveyors with a full length enclosure. In order to decrease the self-heating effect in the silo and prevent fire, two blowers each 40 hp blow air from the fully perforated bottom of the silo.  The temperatures inside of the silo are monitored constantly using 12 cables suspended from the top at different locations with sensors at 10 levels for each cable. The silo dimension and cable location is illustrated in Figure B.1 (Appendix B). 123  Figure 7.1 Illustration of the cylindrical silo and two dimensional axi-symmetric model 124  We chose one silo to simulate the self-heating process. Assuming the system is axi- symmetric, a two dimensional axi-symmetric model is developed in this simulation. Axi- symmetric means variations in temperatures are along the radius and the height of a silo. The factors studied include airflow rate, wall insulation, ambient wind condition, ambient temperature and storage load. 7.2. Computation domain and boundary condition A two dimensional axi-symmetric computational domain is shown in Figure 7.2. The cylindrical storage domain has a diameter (D) and height (H=D). The silo is 80% full of pellets and the wall thickness is set as (Lw=D/100). The properties of wood pellets are assumed to be homogeneous in this model and are specified in the model section. COMSOL 3.5a software is used to solve the model equations for both models with and without ventilation air flow. The boundary conditions of both cases are shown in Figure 7.2. The left boundary (B.1) is set at the axis (r=0) with a symmetrical boundary condition for both cases. At the right boundary (B.3), the heat flux on the outside silo wall is calculated by  , inf ( )wall Wall W SQ h A T T   (7.1) where QWall is the heat flux from the silo wall to the environment, W. hWall is the heat transfer coefficient between the wall and air W/(m2 K). Eq.(7.1) is applied to both the top and the bottom boundaries as well. The interior boundary between the wall and bulk pellets is set by the continuity equation, expressed as  P P W WT T      . (7.2) 125  where the left side indicates the heat flux of the bulk pellets on the boundary while the right side indicates the heat flux across the wall or the boundary. λP is the thermal conductivity of wood pellets, W/(mK); λW is the thermal conductivity of the wall, W/(mK). Case A : no airflow B. 1 B. 3 B. 2 B. 4 Case B : with airflow B. 1 B. 4 B. 3 B. 2 Boundary Condition: B.1 Axial symmetry at r=0. B.2 Newton cooling law Eq.(7.1) at boundary B.2, B.3, B.4 (outside walls) Boundary Condition for Air: B.1 Axial symmetry at r=0. B.2 Airflow inlet, v=v0 ; Tair=T_inf B.3 No slip wall. B.4 Airflow outlet, constant pressure Boundary Condition for Soild: B.1 Axial symmetry at r=0. B.2 Thermal insulation, B.3 Newton cooling law Eq.(7.1) B.4 Thermal insulation, R=D/2 L W H LW  Figure 7.2 Diagram of computation domain and boundary condition for both cases with and without airflow 126  7.3. Governing equations and computational module 7.3.1. Governing equations The energy balance for a cylindrical storage can be describe by the following equation For the solid phase (Wood Pellets)  2 1 1 1 1 ,1 1 12 1 ( ( ) ) ( )P G H T T T C r Q T Q t r r r z               (7.3) where 1 is the bulk density of the pellets, kg/m3; ,1P C  is the specific heat capacity of the bulk wood pellets, J/(kg K); T1 is temperature of the wood pellets, K; 1  is the effective thermal conductivity of wood pellets, W/(m K); QG(T1) is the heat generation from the bulk pellets, W/m3; QH is the convective heat flux from wood pellets to airflow, W/m 3; r and z are the cylindrical coordinates, m; t is the time, s. For gas phase (Air) Energy 2 2 2 2 2 2 ,2 2 2 1 ( ) ( ( ) )P z H T T T T C r Q t z r r r z                  (7.4) Momentum 2 2 2 22 1 ( ) ( ( ) )z z z zz z P r g t z r r r z z                             (7.5) Continuity 2 2 ( ) 0z t z          (7.6) where 2 is the bulk density of the air, kg/m3; ,2P C  is the specific heat capacity of air, J/(kg K); T2 is temperature of air, K; 2  is thermal conductivity of air, W/(m K); z is advection velocity of the gas phase in z direction, m/s; µ2 is the viscosity of air, (Pa s); QH is the convective heat flux from wood pellets to airflow, W/m3; r and z are the cylindrical coordinates, m; t is the time, s. ϕ is the porosity of the pellets, in fraction.   127  7.3.2. COMSOL modules COMSOL is used in this study to solve model equations. Three modules are used in this simulation to solve the coupled energy and momentum equations. Heat transfer by conduction (ht) is used to solve the solid phase heat transfer equation, Eq.(7.3). Convection and conduction (chcc) and Incompressible Navier-Stokes (chns) module are used to solve the gas phase heat and momentum transfer equations[130], Eqs. (7.4-6). Heat Transfer by Conduction (ht) (Solid phase) Convection and Conduction (chcc) (Gas phase) Incompressible Navier-Stokes(chns) (Gas phase) Output {‘T’,’T2',’v’,’p’} Input {‘R’,T_init’,’T_i nf’’,’v0'} dQ=hloc(aSdz)(T-T2) {‘T2’,’v’,’p’}  Figure 7.3 Diagram of the COMSOL modules used in this simulation 7.4. Model parameters In this section, the selection of the model parameters is described with the complete material properties that used in the simulation given in Appendix A. 7.4.1. Wood pellets properties The bulk density of the wood pellets ranges from 600 to 700 kg/m3 and the average of 650 kg/m3 is used in this study. The thermal properties of wood pellets (λ1 and ,1PC ) have been 128  investigated in the previous chapters and the results show that the following equations can be used to estimate the properties of wood pellets.  1 (0.219 0.01 ) (1 ) 0.027M         (7.7)  P,1C = (1- M /100)(5.46T -524) + 42M  (7.8) where M is the moisture content of the pellets, in percentage wet basis; λ1 and CP,1 are the effective thermal conductivity, in W/(mK) and specific heat capacity J/(kg K) of the wood pellets. 7.4.2. Heat generation source In Chapter 3, the heat generation rate of wood pellets was measured experimentally at low temperatures (30oC- 50oC). The results showed that small amount of heat was released during wood pellets storage due to slow oxidation and the heat generation rate was highly dependent on the pellets age and storage temperature, while the pellets moisture content had almost no affect on the heat generation. Eq. (7.9) was developed to estimate the heat generation rate of wood pellets in a sealed environment and will be used in the current simulation.  G 1 1 11637.80 ( ) 1.15 15exp( 0.0256 ) exp[ 2.883exp( 0.014 6801.78 / ) 1 4 )] 0.0026 Q T e x T x T e t               (7.9) Because the oxidation reaction is extremely slow at low temperatures (30oC -50oC), it is reasonable to assume the presence of sufficient oxygen in the voluminous and open silo for the low temperature oxidation and Eq. (7.9)becomes  1G 11637.80 ( ) 1.15 15exp( 0.0256 ) 0.0026Q T e x T        (7.10) where QG is the heat generated from wood pellets, W/m 3. ρ is the bulk density of wood pellets, kg/m3. x1 is the age of the pellets, day, and equal to x0+t/86400, where x0 is the initial age of the pellets, day, and t is time in second. 129  In Chapter 5, the heat generation rate of wood pellets at high temperatures (100oC to 200oC) was experimentally measured. The experimental results were analyzed using both analytical method and numerical method. Eq. (7.11) was obtained for estimating the heat generation rate before ignition.  9( ) 4.22 10 exp[ 78700 / (8.314 )]GQ T T      when 323K<T<473K (7.11) where, QG is the heat generated from wood pellets, W/m 3. ρ is the bulk density of wood pellets, kg/m3. T is temperature, K. 7.4.3. Heat transfer on the outer wall surface The dimensionless heat transfer coefficient (Nusselt number) depends on the property and the velocity of the fluid and can be expressed as functions of Reynolds number and Prandtl number. Considering the general situation, the heat transfer at the outside wall depends on the wind condition, which is less than 1 km/h for calm air and 40 km/h for storm. Both free convection and forced convection are considered in this study. Twall Tinf Tinf Twall Free Convection Force convection Wind condition H D Pellets Pellets  Figure 7.4 Diagram of the airflow and temperature pattern The mean Nusselt numbers (Num) is calculated base on two airflow patterns, as shown in Figure A.1. When there is no wind outside of the silo wall, the free convection creates an 130  vertical airflow from the bottom upward and a thin laminar boundary layer is assumed near the wall as shown in Figure 7.4a. For conditions with wind flow, a horizontal airflow across the silo is assumed as shown in Figure 7.4b. 7.4.3.1. Free convection Free convection is considered during no wind condition and a thin laminar boundary is assumed near the wall (as shown in Figure 7.4a). The Nusselt number can be calculated following [131]  1/4(Pr, )( Pr)freeNu C shape Gr  (7.12) where, C is a function of Pr and the shape of heat transfer surface and can be expressed as  1 2( ) (Pr)C C shape C  (7.13) Gr is the Grashof number. The Prandtl number of the air is calculated to be 0.70-0.77 for a temperature range of 273-473 K. C1 was found to be 1.0 for vertical plate and 0.835 for horizontal plate (top boundary); C2 was found to be 0.515 [131]. The second term on the right side of Eq. (7.12) can be written as  2 3 0 1 2 ( ) ( Pr) ( )( )P g H T T C Gr         (7.14)  where,  is the gas density at temperature T . β is the coefficient of volume expansion which is defined as 1 1 ( ) ( )P P V V T T           . Inserting the air properties, the heat transfer coefficient can be expressed as  2 3 1/4 1/4 ,1 1 2 3 4[( ) ] [( ) ] P Wall free gCk T h Nu C C C C H TH        (7.15) 131  where C3 is a function of  air properties; C4 is a function of H (height of the silo) and ∆T (temperature difference between wall and ambient air). H is replaced by D (silo diameter) for heat transfer on the top surface of the silo. Figure 7.5 shows the estimated free convection heat transfer coefficient at different wall temperatures and silo heights. At a constant ambient temperature (288 K), the heat transfer coefficient is significantly affected by the silo height and the temperature difference between wall surface and ambient tempreatures. When the wall temperature is 303 K, the heat transfer coefficient of 1 meter silo is 2.7 W/(m2K), and that of 20 meter silo is 1.3 W/(m2K). The heat transfer coefficient is higher when the temperature difference between the wall temperature and ambient temperature is larger. At an ambient temperature of 288 K, the heat transfer coefficient of 1 meter silo increases from 2.2 to 3.9 W/(m2K) when the wall surface temperature increases from 295 K to 373 K. 290 300 310 320 330 340 350 360 370 380 0.5 1 1.5 2 2.5 3 3.5 4 Wall surface temperature (K) H e a t tr a n s fe r c o e ff ic ie n ts  ( W /( m 2  K ))   H=1m H=5m H=10m H=20m  Figure 7.5 Estimated heat transfer coeffient between the wall surface and the amibent air (temperature 288K) under free convection condition. Note: H is the silo height .  7.4.3.2.  Force convection (Wind condition) In addition to the vertical airflow caused by the natural convection, a force convection caused by wind is also considered (as shown in Figure 7.4b). The wind direction is assumed 132  to be horizontal. When the wind speed ranges from1 km/h to 40 km/h (0.2-11m/s), the Reynolds number is in the range of 3.4×105-136×105 for a 20 meter diameter silo. In this flow range, Nu can be calculated using the following empirical correlation  1/2 1/30.6Re Prconv freeNu Nu   (7.16) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 1 2 3 4 5 6 7 8  Wind velocity outside of silo (m/s) H e a t tr a n s fe r c o e ff ic ie n ts  ( W /( m 2  K ))   D=1m D=5m D=10m D=20m  Figure 7.6 Estimated heat transfer coeffcient between the wall surface (305K) and the ambient air (295K) under wind condition. Note: D is the diameter of the silo Figure 7.6 shows the estimated heat transfer coefficient between the wall surface (305K) and the ambient air (295K) at different wind conditions. The heat transfer coefficient is dependent on the silo diameter and the wind velocity. For a 1 m diameter silo, it increases from 4.0 W/(m2K) to 7.5 W/(m2K) when the wind speed increases from 0 to 2 m/s. The heat transfer coefficient for 20 m silo is much smaller and increases from 1.2 W/(m2K) to 2.4 W/(m2K) when the wind speed increases from 0 to 2 m/s. The wind effect on the heat transfer coefficient is smaller for larger silos. As the characteristic length of the silo increases, the dimensionless heat transfer coefficient (Nu) increases; however, the heat transfer coefficient (h) decreases. The heat transfer coefficient of gas can also be found in many literatures, 3-20 W/(m2K) for free convection and 10-100 W/(m2K) for forced convection [112, 131]. The heat transfer coefficient estimated from the empirical equations seems to be smaller than the literature values because 133  of the large characteristic length of the silo. A sensitivity analysis is included in the Appendix C to discuss the effect of the heat transfer coefficient on the simulation results. 7.4.4. Heat exchange between airflow and wood pellets inside of silo Heat transfer coefficient hloc between wood pellets particles and air flow in the packed bed is an important parameter in the heat transfer model. The heat transfer coefficient is highly dependent on the velocity profiles of the air flow and the void fraction of the packed bed. hloc is defined for a representative volume Sdz of the packed bed following [131] d z Wood Pellets (T1) Air flow (T2) dQ   1 2( )( )H loc SdQ h A Sdz T T   (7.17) where QH is the convective heat flux from wood pellets to airflow, W/m 3. AS is the outer surface area of wood pellet particles per unit bed volume or specific surface area, m2/m3(packed bed), and hloc is the local heat transfer coefficient between pellets and airflow, W/m2 K. The forced convection for gas flow through packed beds has been analyzed widely and the local heat transfer flux correlation has been reported as  2/3 0.3812.19Re 0.78ReHj     (7.18) Where jH, the Chilton-Colburn factor, and Reynolds number are defined as  2/3 0 ( )loc PH P h C j C G k    (7.19) 134   0 0 6 Re (1 ) S S D G G A       (7.20) Here, 0G   is the superficial mass flux; ν is the inside airflow speed. The quantity ψ is the particle-shape factor, with a fitted value of 0.92 for cylindrical pellets. AS is the outer surface area of the particles per unit bed volume, m2/m3(packed bed). Ds is the particle diameter, m and ϕ is the porosity in fraction. Elimination of jH and Re using the three equations (Eq. (7.18), (7.19) and (7.20))gives  1/3 2/3 0.619 0.381 0.286 2 1/3 0 00.663( ) ( ) 0.394 ( ) ( )loc P S S Ph G C A k G A C k      (7.21) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 20 40 60 80 100 120 140 Velocity of the forced ventilation inside silo (m/s) C o n v e c ti v e  h e a t tr a n s fe r c o e ff ic ie n ts  ( W /( m 2 K ))   T=273K T=373K T=473K  Figure 7.7 Estimated convective heat transfer coefficient vs. gas ventilation velocity at different temperatures The convective heat transfer coefficients at airflow velocity of 0.001-1 m/s were estimated at temperatures of 273K, 373K and 473K as shown in Figure 7.7. The heat transfer coefficient is highly dependent on the airflow speed and less on the temperature; it was estimated to be 4.8 to 123 W/(m2K) at 273 K when the airflow increased from 0.001-1 m/s. Eq. (7.21) was used in the model to estimate the convective heat transfer coefficient between the airflow inside the silo and the wood pellets particles. 135   7.4.5. Air flow velocity generated by the fans Air injection from the bottom of the silo is one of the most effective methods to reduce the self-heating effect and to cool the wood pellets down. The air velocity of the flow is calculated in order to make the simulation close the real case. For the Fiberco silos, there are two fans of 40 horsepower at the bottom of each silo. By the definition of the mechanical power  W v p S  (7.22) where W is the power of the fans, 59680 W (40*2 hp); ν is the velocity of the airflow, in m/s ; ∆p is the pressure drop along the silo height, in Pa; S is the cross sectional area of the cylinder, in m2. Yazdapanah [29, 30] developed the relationship between the airflow velocity and pressure drop based on measurements. The experimental result showed that the pressure drop followed the Ergun equation,  2p a b L       (7.23) where 2 3 2 2 (1 ) 150a d       2 2 1 1.75b d      where ϕ is the porosity of the packed bed, in fraction; ρ is the bulk density, in kg/m3; d is equivalent diameter of the pellet particles, in m; μ is the viscosity of the air, in Pa∙s; and ψ is the shape factor of the pellets. 136  The airflow superficial velocity ν (m/s) can be solved by combining Eq. (7.22) and Eq. (7.23). The results are shown in Figure 7.8. The superficial airflow velocity is a function of the packed height of the pellets inside of the silo and the fan efficiency. The airflow speed is smaller for a larger load (height) of pellets inside the silo. Assuming a fan efficiency of 80% the airflow speed ranges from 0.21 to 0.049 m/s for pellets (load) height of 1 to 19 m inside of the silo. For a 80% full silo (14 m height), the superficial airflow speed in the silo is 0.59 m/s when the fan efficiency is 80%. 0 2 4 6 8 10 12 14 16 18 20 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22  Height of the bulk pellets (m)  I n s id e  a ir fl o w  v e lo c it y  ( m /s )   20% effciency 40% effciency 60% effciency 80% effciency  Figure 7.8 Estimated airflow velocity inside of the silo (Fiberco, Inc.) 7.4.6. Mesh size and time step A silo of 20 meter in diameter (H=14 meter, LW= 0.1 meter) is used in this simulation to investigate the effect of mesh size and the time step. An average ambient air flow velocity of 10 km/h (2.8 m/s) is considered in this study (Re=1.8×106) and Eq. (7.16) is used to estimate the heat transfer coefficient h. 137  60s 360s 3600s 2h 4h 6h 12h 0 50 100 150 200  D if fe re n c e  Timestep  a.First stage, T init=300K 60s 360s 3600s 2h 4h 6h 12h 24h 0 0.2 0.4 0.6 0.8  D if fe re n c e  Timestep  b.Second stage, T init=300K 60s 360s 3600s 2h 4h 6h 12h 24h 0 10 20 30  D if fe re n c e  Timestep  c.Third+Fourth stage, T init=300K  Figure 7.9 The residues of different stages In order to select the proper time step for the simulation, the time steps of 60 s to 24 h were tested in each of the stage (indicated in 6.5). Assuming that the smallest time step gives the most accurate results (T0), the lager time step will decrease the accuracy of the results. The residue of each time step is defined by the difference between the results of each time step and the smallest time step size (T0), and expressed as  2 1 0 1 0 0 ( ) t end r R z H t r z Residue T T          (7.24) 138  The residues of different stages are listed in Figure 7.9. The residues start to increase visibly, especially at the first stage, therefore a time step of 3600s is used in the calculation. A triangular mesh automatically generated from COMSOL was used and the mesh size is also tested and the maximum element size used is 1/20th of the radius. 7.4.7. Criterion of thermal runaway Unlike in the experiment where the thermal runaway (pyrolysis, smoldering or ignition) can be observed visually, a criterion need to be developed to define when the runaway happens in the mathematical simulation. As we know from the previous chapter, pyrolysis becomes significant and temperature rises rapidly when the temperature inside the pellet pile exceeds 200 oC. Also, the critical ambient temperature for the small container is 191 oC. Therefore, in the mathematical model, a value of 200 oC was selected as the criterion to indicate a thermal runaway. 7.5. Evolution of self-heating and thermal runaway Figure 7.10 shows the evolution of self-heating of three cases, ambient temperature of 310 K, 315 K and 320 K for a silo of 10 m in diameter and the ambient wind speed of 1 m/s. According to the model results, self-heating is developing very quickly at the beginning and becomes slower after 10 days. In the first stage, the slow oxidation reaction first happens on the low molecular weight substances which are most susceptible to oxidation. When those substances run out (point Ae), self-heating is developing slightly slower in the second stage. When the temperature reaches around 100 oC (point Be), pyrolysis start to take place and the self-heating becomes dangerously rapid (Ce). The heating rate in this stage is so large that it quickly leads to ignition in the center of the silo (over 650 K) where the temperature is seen to shoot up (De). Therefore, we define the third stage (Be →Ce) as the thermal runaway zone, where the self-heating rate increases rapidly for the second time. 139   0 10 20 30 40 50 60 70 80 300 320 340 360 380 400 420 440 460 480 500 Time (days) C e n te r te m p e ra tu re  ( K )   T  init=310K T  init=315K T  init=320K   r,meter z ,m e te r A e  Time=3days   0 5 0 2 4 6 300 350 400 450 500 r,meter z ,m e te r B e  Time=10days   0 5 0 2 4 6 300 350 400 450 500 r,meter z ,m e te r C e  Time=20days   0 5 0 2 4 6 300 350 400 450 500 r,meter z ,m e te r D e  Time=25.5days   0 5 0 2 4 6 300 350 400 450 500  Figure 7.10 Four stages of self-heating. 140  For the case of an ambient temperature of 320 K, the temperature contours at the end of each stage are shown in Figure 7.10. At the end of first stage (Ae), the maximum temperature reached 335 K and the temperature distributes quiet evenly in the silo with a small gradient near the wall. This means that low temperature self-heating is developed more uniformly in the silo due to the slowness of the self-heating reaction at low temperatures. At the end of stage two (Be), the temperature rises to 360 K in 20 days with a larger temperature gradient in the silo. At the end of stage three (Ce), the temperature reaches 380 K in 10 days and the temperature gradient grows larger. At stage four (De), the temperature finally shoots up with an extremely high temperature gradient shown in the graph. This large temperature gradient in stages three and four indicates that at these two stages, the heat generated from the pyrolysis reaction is too large to be dissipated, leading to the accumulation of heat in the silo center. 7.6. Parametric study Considering that the wood pellets are stored in a cylindrical silo, the parameters of the storage condition are studied in this section. 7.6.1. Thickness of the wall To investigate the wall insulation effect on the heat dissipation, the material and the thickness of the wall are examined during the self-heating process. The thickness of 1 cm to 15 cm for steel and concrete were examined in this study. The wall thickness has almost no effect on the temperature at the axis and 8 m away from the centre, but has some effect on the temperature on the wall and in the near wall region, as shown in Figure 7.11. Figure 7.11a and b show that the outer wall temperature and inner wall temperature for concrete and steel wall. The thickness of the wall have negligible effect on the wall temperature because of the good conduction of the steel material. For concrete wall, the outer wall temperatures increase with the wall thickness while the inner wall temperature remains the same. This is caused by the smaller thermal conductivity of the concrete materials which causes the slower heat dissipation. 141   0 10 20 30 40 50 60 70 300 301 302 303 304 305 Time, days T e m p e ra tu re ,K a. Outside wall temperature of different thickness and material   0 10 20 30 40 50 60 70 300 300.5 301 301.5 302 302.5 303 Time, days T e m p e ra tu re ,K b. Inside wall temperature of different thickness and material   Lc=0.025m, Concrete Lc=0.025m, Steel Lc=0.1m, Concrete Lc=0.1m, Steel Lc=0.15m, Concrete Lc=0.15m, Steel Lc=0.025m, Concrete Lc=0.025m, Steel Lc=0.1m, Concrete Lc=0.1m, Steel Lc=0.15m, Concrete Lc=0.15m, Steel  Figure 7.11 The wall thickness and wall material effects on the outside and inside wall temperature.  7.6.2. Wind speed outside of the wall The wind speed of 0-3 m/s is used in the simulation for steel wall silos (20 meter diameter) at ambient temperature 310K. The results are given in Figure 7.12 and Figure 7.13. The temperatures at three points, center (r=0 m, z=7 m), 8 m (r=8 m, z=7 m), 9 m (r=9 m, z=7 m), are presented in Figure 7.12. The results show that the center and 8 m location temperatures are almost identical at different wind speeds (0, 1, 2, 3 m/s), while there exists a small difference at the 8 m location. In another word, the wind speed only affects the 142  temperature really close to the wall, but not on the center regions. In order to further investigate the wind effect in the region near the wall, the inner wall surface temperature and outer wall surface temperature are plotted in Figure 7.13. The results show that the fast the wind speed is; the lower the wall surface temperature is for both inner and outer surfaces. 0 10 20 30 40 50 60 70 310 315 320 325 330 335 340 345 350 355 Time (days) T e m p e ra tu re  ( K )    r=0m  r=8m  r=9m  Figure 7.12 The temperature profile of steel wall at different wind speeds v=(0,1,2,3 ) m/s 0 10 20 30 40 50 60 70 309.5 310 310.5 311 311.5 312 312.5 313 313.5 314 Time (days) T e m p e ra tu re  ( K )    Inside wall  Outside wall v=0 m/s v= 1 m/s v= 2 m/s v= 3 m/s  Figure 7.13 The temperature development of inner and outer walls atdifferent wind speeds 143  7.6.3. Dimension of the container As discussed in Chapter 6, the dimension of the container has a significant effect on the thermal runaway and critical temperature for spontaneous combustion. The temperature development in containers of radius 0.2-10 meters is predicted and the critical ambient temperatures corresponding to thermal runaway are shown in Figure 7.14. Comparing with Figure 6.8, the predicted critical ambient temperature for thermal runaway is much lower in Figure 7.14. This lower temperature is likely associated with the released heat from self- heating at low temperatures. In another word, the heat generation at low temperatures strongly decreases the critical ambient temperature. When the ambient temperature is higher than the critical ambient temperature (in Figure 7.14), a thermal runaway will possibly happen. Therefore, the area under the critical ambient temperature line can be considered as a safe region. The critical ambient temperature decreases when the radius of the silo increases. The critical ambient temperature for a 10 m radius silo is predicted to be 309 K, or 36 oC. 0 2 4 6 8 10 300 320 340 360 380 C ri ti c a l a m b ie n t te m p e ra tu re  ( K ) Radius (m) Spontaneous Combustion Safe Region  Figure 7.14 Relationship between critical ambient temperature and radius of the silo 144  Note: Critical ambient temperature is defined as the minimum ambient temperature that can lead to a thermal runaway. 7.6.4. Heat transfer model with forced air flow In non-isothermal flow, the fluid properties change with the fluid temperature and pressure. The momentum change is equal to the pressure change in addition to the forces along the flow direction.  The momentum equation, Eq. (7.5), including both force convection and natural convection, can be expressed as[112] [131].  2 2 2 2 22 1 ( ) ( ( ) ) ( )z z z zz z P r g g T T t z r r r z z                                 (7.25) where 2 is the bulk density of the air, kg/m 3; υz is advection velocity of the gas phase in z direction, m/s; µ2 is the viscosity of air, (Pa s); r and z are the cylindrical coordinates, m; t is the time, s. ϕ is the porosity of the pellets, in fraction.   is the average thermal expansion coefficient, 1/K. The momentum changes during forced convection are primarily determined by the external force ( 2 z P g z      ) and buoyancy force ( 2 ( )g T T   ). The buoyancy force caused by the temperature gradients, is relatively small [131] as compared to the external convective force. Therefore, it is safe to ignore the buoyancy force during forced convection and Eq. (7.25) can be simplified to Eq. (7.5). However, if there is no forced convection the buoyancy force ( 2 ( )g T T   ) becomes the primary force that could affect the flow pattern. Figure 7.10 shows that the temperature profiles inside the wood pellets silo developed very slowly and evenly before the thermal runaway. The airflow caused by the buoyancy force, 2 ( )g T T    was less than 10 -6 m/s, which was negligible as compared to the forced air velocity (0.0001-1m/s). Therefore, in this section, which focused on the effect of the forced air flow, the free convection effect was neglected. The effect of free convection without forced air flow is discussed in Appendix E. 145  7.6.4.1. Continuous airflow Assume vertical air flow evenly distributed from the bottom, the airflow effect is simulated using boundary condition for Case B as described in Figure 7.2 and COMSOL Multiphysics modules described in Figure 7.3. The airflow velocity range of 0.0001-0.3 m/s is investigated in a silo diameter of 21 m and a height of 17 m (80% full or loaded). From Figure 7.14, we know that the critical ambient temperature is 309 K for this silo. When the ambient temperature is higher than 309 K, a thermal runaway is likely to occur in the silo. Figure 7.15 shows some simulation results of the airflow cooling effect at several ambient temperatures. The airflow cooling effect is significant when the ambient temperature is low (314 K and 319 K). At an ambient temperature of 314 K, airflow cooling effect at the air velocity range of 0.0001-0.02 m/s is illustrated in Figure 7.15a. A thermal runaway happens when the air velocity is set at 0.0001m/s. When the airflow is higher than 0.001 m/s, the temperature at the centre goes up at the first 10 days but eventually decreases and no thermal runaway is detected. Similar trend is observed in Figure 7.15 b, where a thermal runaway is observed when airflow is at 0.001m/s but no thermal runaway is observed for airflow set as 0.002 m/s. At a higher ambient temperature (329 K, 335 K), the airflow effect is limited and a higher velocity is required to prevent the thermal runaway as illustrated in Figure 7.15c,d. An airflow velocity of 0.055 m/s is required for an ambient temperature of 329 K and an airflow velocity of 0.24 m/s is required for an ambient temperature of 335 K.  146  0 10 20 30 40 50 60 70 80 310 320 330 340 350 360 370 380 390 400  T e m p e ra tu re  ( K )  Time (days) results 314K   v=0.0001 m/s v=0.001 m/s v=0.002 m/s 0 10 20 30 40 50 60 70 80 310 320 330 340 350 360 370 380 390 400  T e m p e ra tu re  ( K )  Time (days) results 319K   v=0.001 m/s v=0.003 m/s v=0.004 m/s  a. Ambient temperature 314 K b. Ambient temperature 319 K 0 2 4 6 8 10 320 330 340 350 360 370 380 390 400 410 420 430 results 329K  T e m p e ra tu re  ( K )  Time (days)   v=0.02 m/s v=0.04 m/s v=0.055 m/s 0 2 4 6 8 10 320 330 340 350 360 370 380 390 400 410 420 430 results 335K  T e m p e ra tu re  ( K )  Time (days)   v=0.19 m/s v=0.24 m/s v=0.29 m/s  c. Ambient temperature 329 K d. Ambient temperature 335 K  Figure 7.15 Simulation results of airflow effect at different ambient temperatures and airflow velocities A3 A4 A1 A2 147  r (m) z  ( m ) A1 T=314k, v=0.0001m/s Time=10days   0 5 10 0 5 10 320 340 360 380 400 420 r (m) z  ( m ) A2 T=314k, v=0.0001m/s Time=62days   0 5 10 0 5 10 320 340 360 380 400 420 r (m) z  ( m ) A3 T=329k, v=0.02m/s Time=2days   0 5 10 0 5 10 320 340 360 380 400 420 r (m) z  ( m ) A4 T=329k, v=0.055m/s Time=3days   0 5 10 0 5 10 320 340 360 380 400 420  Figure 7.16 Solid phase temperature profiles at different airflow rates Figure 7.16 illustrates the simulated solid phase temperature profiles at different days with different airflow rates. For instance, Figure 7.14a shows the solid phase temperature profile at 10 days at an ambient temperature of 314 K and an airflow velocity of 0.0001 m/s. Because the velocity is really slow, the temperature at 10 days (in Figure 7.14a) is very similar with one without airflow (Figure 7.10) where the temperature gradients are mainly near the wall. Figure 7.14b shows the temperature profile while thermal runaway (at 62 days). While the highest temperature in Figure 7.10 (temperature profiles without airflow) was at the center of the silo, the highest temperature in Figure 7.16b (temperature profiles with airflow) was at the top of the silo. The temperature vertical gradients inside of the silo increase with the airflow velocity as shown in Figure 7.16c and Figure 7.16d. Figure 7.16c and Figure 7.16d show the simulation results of point A3 and A4, which are indicated in Figure 7.15. This is because the airflow exchange heat with pellets along the path inside the silo so that the heat is removed from the top by the airflow. This may explain why fire often started from the top of the silo in some cases. 148  As we observed, the higher the ambient temperature is, the higher air velocity is required to prevent the thermal runaway inside the silo. The minimum velocity which can prevent the thermal runaway is defined as the critical airflow velocity. The critical airflow velocity at different ambient temperatures is calculated for the silo of 21 m in diameter and shown in Figure 7.17. The airflow is seen to have a significant cooling effect when the ambient temperature is lower than 330 K and the critical airflow velocity is found to be 0.05 m/s. Therefore, a continuous airflow of 0.1m/s as used by Fiberco Inc. is sufficient to prevent thermal runaway for ambient temperatures lower than 330 K (57 oC). However, the airflow cooling effect is limited if the ambient temperature is higher than 330 K (57 oC). 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 310 320 330 340 350 C ri ti c a l a m b ie n t te m p e ra tu re  ( K ) Airflow velocity (m/s) Safe Region Spontaneous Combustion  Figure 7.17 Relationship between ambient temperature and airflow velocity 7.6.4.2. Periodical airflow Besides continuous airflow cooling, periodical airflow ventilation is also investigated in this study. Base on the same model equations (boundary condition for case B and COMSOL modules described in Figure 7.3), a periodical airflow cooling effect is studied. The airflow is issued every other 24 hours. The results are shown in Figure 7.18. 149  0 10 20 30 40 50 60 310 320 330 340 350 360 370 380 390 400 Time (days) C e n te r T e m p e ra tu re  ( K )   v=0.001 m/s v=0.002 m/s v=0.003 m/s  0 5 10 15 20 25 30 310 320 330 340 350 360 370 380 Time (days) C e n te r T e m p e ra tu re  ( K )   v=0.004 m/s v=0.005 m/s v=0.006 m/s  a. Ambient temperature 314K, airflow every other 24 hours b. Ambient temperature 319K, airflow every other 24 hours Figure 7.18 Solid phase temperature development with periodical airflow Comparing with Figure 7.15 which shows the continuous airflow, periodical airflow is less efficient than the continuous airflow. For example, at an ambient temperature of 314 K (41 oC), a continuous airflow of 0.001 m/s is required to prevent the thermal runaway (Figure 7.15) while a periodical airflow of 0.002 m/s is required in Figure 7.18a. At an ambient temperature of 319 K, a continuous airflow of 0.003 m/s is required to prevent the thermal runaway (Figure 7.15b) while a periodical airflow of 0.005 m/s will not be able to prevent the thermal runaway Figure 7.18b. From the results of airflow simulation (Figure 7.15 and Figure 7.18), it is found that the thermal runaway normally occurs in the first 10-20 days of storage because the self-heating reaction activity become lower when the pellets getting old according to Chapter 4 and Chapter 5. After 20 days, self-heating from low temperature oxidization reactions will phase out if the temperature is not high enough to move to the next stage.  150  7.7. Model verification 7.7.1. Comparison with industrial silos For the purpose of verification, temperature data were collected from a real industrial silo as was explained earlier in this chapter. Each silo was equipped with 12 cables and each cable had 11 sensors located at different levels of the silo as shown in Figure B.1 (Appendix B). Assuming an axi-symmetric silo, Figure 7.19 shows the temperature profile of cable 1 (11 levels). In this data set, there was no airflow in the first 5 days and the airflow started from the day six. During the first five days, the center level temperature (L6 and L7) rose up the fastest and reached to 44 oC on the fifth day. The hot spot located at the center at no airflow, consistent with what predicted in Figure 7.10. After day six, the airflow started to be issued into the silo from the bottom and the hot spot started to move up. The temperature at the top (L11) developed to the highest temperature (57 oC) inside the silo and the peak temperature at each level decreased as their level decreased. The bottom level (L1) had the lowest temperature. This phenomenon is in agreement with the temperature patterns predicted in the simulation results (Figure 7.16). The fan was controlled by the relative humidity of the air, which was automatically turned off when the relative humidity was over 90%. In the simulation, we used a periodic airflow of every other 24 hours after the sixth day. The simulation results for air velocities of 0.06 m/s and 0.03 m/s are shown in Figure 7.20. 151  0 5 10 15 20 25 10 15 20 25 30 35 40 45 50 55 60  L11  L10  L9  L8  L7  L6  L5  L4  L3  L2  L1T e m p e ra tu re , o C days from 10th June, 2009  Figure 7.19 Temperatures of cable C1 inside a industrial silo (80% load) 0 5 10 15 20 25 30 310 315 320 325 330 335 Time (days) C e n te r te m p e ra tu re (K )    13meter  8.6meter  6.4meter  4.2 meter  Ground 0 5 10 15 20 25 30 310 315 320 325 330 335 Time (days) C e n te r te m p e ra tu re (K )    13meter  8.6meter  6.4meter  4.2 meter  Ground  a.V=0.06 m/s b.V=0.03 m/s Figure 7.20 Simulation results for airflow velocities of 0.06 m/s and 0.03 m/s. 152  Figure 7.20 shows the self-heating simulation results of the commercial silo at ambient temperature of 37oC and inside airflow velocity 0.06 m/s and 0.03 m/s. In Figure 7.20, the axial temperature profiles at 5 axial levels are plotted, 13 m, 8.6 m, 6.4 m, 4.2 m (above the ground) and ground. At the initial 5 days, the temperatures at all levels reach 326 K (53oC). After the sixth day, the airflow begins to be issued to remove the heat away from the top opening. When the air velocity is set at 0.06 m/s, the temperatures decrease quickly and become stable after 5 days of ventilation. The peak temperature reaches 330 K (57oC) at the 7th day. When the air velocity is set at 0.03 m/s, the temperatures decrease much slower and stabilize after 5 days of ventilation. The peak temperature reaches 332 K (59 oC) with the top being the highest. Because the silo is 80% full, top surface of the packed solids (pellets) is at around level 8. Comparing with the industrial data Figure 7.19, the temperature profiles of Figure 7.20a are very close to the industrial data. We assumed an evenly distributed airflow and a two dimensional axi-symmetric system in the simulation. However, the real airflow is not evenly distributed, which will cause a systems error to the simulation results. 7.7.2. Pilot scale experiments A pilot experiment has been conducted in a cylindrical reactor of 1.2 m in diameter and 4.6 m in height, as shown in Appendix C (Figure C.1). The reactor is equipped with 30 thermocouples located at 6 different levels. After loaded 3 metric tons fresh wood pellets, the reactor was sealed and temperatures were monitored for 60 days. Figure 7.21 shows the temperature profiles at level 1 (lower level) and level 3 (medium level). At the lower level (Figure 7.21a), the initial temperatures of the bulk wood pellets was 18.5oC and the ambient temperature (room temperature) was about 21oC; both center temperature and the half way temperature (half of the radius to the center) increased with time and reached the wall temperature after about 80 hours. After reached the wall temperature, they kept rising and reached a maximum temperature about 2oC higher than the wall temperature. Even though the temperature rise is small, the center temperature is slightly higher than the half way temperature (level 1) at the peak, which is an evidence of self-heating in the container. After the peak, the temperature gradually decreased to 22oC at 450 hours, 1oC higher than the wall temperature.. Figure 7.23a shows the simulation results at level 1 with assumed unlimited oxygen supply. The simulation results showed an excellent agreement with the experimental 153  results in the first 160 hours. Both the center temperature and half way temperature increased with time; reached the wall temperature at about 80 hour; kept increasing afterward. However, there is a slight deviation of the simulation from the experiment after 160 hours; the predicted temperature kept rising for another 10 hours instead of going down as observed in the experiment and reached a maximum temperature 0.5oC higher than the experimental results. This deviation might have been caused by the oxygen limitation at the lower level of the silo. The O2 concentration inside the silo was monitored and shown in Figure C.4 (Appendix C). The O2 concentration of the lower level in the silo dropped very quickly with time and reached a value lower than 1% at 160 hours (6.6 days). As a result of self-heating (experiments), the temperature profiles decreased after 160 hours in the experimental results and the unlimited oxygen assumption is no longer valid beyond 160 hours. Figure 7.23a. shows the simulation results using the limited oxygen assumption. The temperatures have the similar trend as in Figure 7.23a (unlimited oxygen supply) but with a lower temperature rise, 1oC higher than the wall temperature, which is the same as the temperature recorded in the experiment at 450 hours. 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level1   Center T7 Half T6 Wall T5 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level3   Center T17 Half T16 Wall T15  a. Temperature profile of level 1 (lower) b. Temperature profiles of level 3 (upper) Figure 7.21 Temperature profiles from pilot scale reactor experimental results. Note: pilot scale reactor, a reactor of 1.2 meter in diameter and 4.6 meter in height (75% full). 154  0 100 200 300 400 500 18 20 22 24 26 28 30  T e m p e ra tu re  (o C )  Time (hours)   Center Half Wall 0 100 200 300 400 500 18 20 22 24 26 28 30  T e m p e ra tu re  (o C )  Time (hours)   Center Half Wall  a. Temperature profile of level 1 (lower) Tambient = 21oC; Tinitial =18.5oC b. Temperature profiles of level 3 (upper) Tambient = 22oC; Tinitial =20.5oC, Figure 7.22 Temperature profiles from pilot size reactor simulation results with unlimited oxygen supply 0 100 200 300 400 500 18 20 22 24 26 28 30  T e m p e ra tu re  (o C )  Time (hours)   Center Half Wall 0 100 200 300 400 500 18 20 22 24 26 28 30  T e m p e ra tu re  (o C )  Time (hours)   Center Half Wall  a. Temperature profile of level 1 (lower) Tambient = 21oC; Tinitial =18.5oC b. Temperature profiles of level 3 (upper) Tambient = 22oC; Tinitial =20.5oC, Figure 7.23 Temperature profiles from pilot size reactor simulation results with limited oxygen supply Ambient temperature was set at 23oC, initial temperature was set at 18.5oC, No ventilation inside or outside. At the upper reactor level (Figure 7.21b and Figure 7.22b), the experimental data (Figure 7.21b) agree well with the simulation results with the assumption of unlimited oxygen supply (Figure 7.22b). Comparing with Figure 7.21a, the temperature profiles at the upper level reached the peak at 200 hours and stabilized afterward. This agrees well with the simulation with the assumption of unlimited oxygen supply (Figure 7.22b). The maximum temperature reached 25oC, which is 3oC higher than the wall temperature. In Figure C.4 (Appendix C), 155  the O2 concentration drop at the upper level (3 m) is slower than at the lower level, reaching 1% over 200 hours (8.3 days). At the upper level close to the top surface (25% head space), there is more oxygen supply than at the lower level. It is reasonable to assume that at the upper level O2 is enough for self-heating over this period. As a result, the simulation results with the unlimited oxygen supply assumption (Figure 7.22b) have a good agreement with the experimental results (Figure 7.21b). More experimental data are given in Appendix C. Due to the natural ventilation in the room, the ambient temperature was generally lower at the lower level of the silo where the fresh air came from the door located near the bottom of the reactor. Therefore, the wall temperature was not uniform and the top air temperature (22oC) was higher than the bottom air temperature (21oC). 7.7.3. Lab scale experiments To verify the simulation at higher temperatures, the thermal runaway in the cylindrical reactor of 417 mm in height and 356 mm in diameter is simulated, with the results shown in Figure 7.24. The reactor was described in Chapter 2. Temperatures were monitored at four locations, T1, T2, T3 and T4, which were 116 mm, 58 mm, 0 mm (center), and 88 mm away from the center, respectively. The experimental results at different oven temperatures were listed in Figure 6.1 and discussed in Section 6.2. In this section, we choose two locations; T1 (116 mm) and T3 (center), to compare between the experimental results and the simulation results. As shown in Figure 7.24, there is a good agreement between the experiment and simulation when the oven temperature was at 90oC (Figure 7.24a) and 105oC (Figure 7.24b). When the oven temperature was 125oC(Figure 7.24c) and 132oC(Figure 7.24d the simulation results showed a slightly higher temperature profile (higher self-heating) than the experiments but no thermal runaway was observed from both experiment and simulation. When the oven temperature was set at 143oC (Figure 7.24e), thermal runaway happened in both simulation and experimental results (spontaneous combustion). Overall, the simulation results agreed well with the experimental data. 156  0 200 400 600 800 1000 1200 1400 1600 20 30 40 50 60 70 80 90 100 Time (min) T e m p e ra tu re  ( o C ) SI11-18-2.xls   Toven T1 experimental results T3 (center) experimental results T1 simulation results T3 (center) simulation results  a. Oven temperature 90oC 0 500 1000 1500 75 80 85 90 95 100 105 110 Time (min) T e m p e ra tu re  ( o C ) SI11-19-2.xls   Toven T1 experimental results T3 (center) experimental results T1 simulation results T3 (center) simulation results  b. Oven temperature 105oC 0 500 1000 1500 2000 2500 3000 95 100 105 110 115 120 125 Time (min) T e m p e ra tu re  ( o C ) SI11-20-2.xls   Toven T1 experimental results T3 (center) experimental results T1 simulation results T3 (center) simulation results  c. Oven temperature 125oC 0 500 1000 1500 2000 2500 3000 120 122 124 126 128 130 132 134 136 138 Time (min) T e m p e ra tu re  ( o C ) SI11-21-2.xls   Toven T1 experimental results T3 (center) experimental results T1 simulation results T3 (center) simulation results  d. Oven temperature 132oC 0 500 1000 1500 2000 2500 3000 3500 120 140 160 180 200 220 240 260 Time (min) T e m p e ra tu re  ( o C ) SI11-24-2.xls   Toven T1 experimental results T3 (center) experimental results T1 simulation results T3 (center) simulation results  e. Oven temperature 143oC  Figure 7.24 Comparison between experimental results and simulation results in lab scale reactor  157  7.8. Conclusions A two dimensional axi-symmetric model was built to simulate the self-heating process and to predict the thermal runaway (spontaneous combustion) during wood pellet long term storage. The self-heating parameters determined in previous chapters were used in this model. A parametric study was carried out to investigate the effects of airflow rate, wall insulation, wind speed, ambient temperature and storage silo dimension. The following conclusions are derived. 1.  Both wall thickness and the convection outside of the storage wall have almost no effect on temperatures of stored pellets, especially the center temperature. They only affect the wall temperature and locations near the wall. Therefore, both wall thickness and convection on the outer wall have a minimal effect on thermal runaway. 2.  Dimensions of the container have a significant impact on thermal runaway and spontaneous combustion. The critical ambient temperatures for silos with a radius of 0.2-10 m are calculated. The results showed that the critical ambient temperature for a 10 m radius silo can be as low as 36 oC. 3.  The inner air ventilation effect is investigated for both continuous airflow and periodical airflow. A critical air velocity is calculated for each given ambient temperature. The results showed that the air ventilation is a very effective approach for reducing self-heating and preventing thermal runaway at ambient temperatures lower than 330 K. The ventilation effect becomes very limited when the ambient temperature is higher than 330 K. 4.  The temperature profiles data from an industrial silo were collected to compare with the simulation results. The simulated temperature developing process from the model is very close to the industrial data. However, the airflow cooling process is faster and more effective in the simulation than observed in the industrial data.    158  8. Conclusions and recommendations 8.1. Conclusions Self-heating is a serious problem during wood pellet storage. This research focused on investigating parameters (conditions) affecting the self-heating process and predicting the critical conditions for spontaneous combustion in stored wood pellets. Experiments were carried out to determine thermal conductivity, specific heat capacity, exothermic properties of wood pellets including heat release rate at both low and high temperatures. A two- dimensional axi-symmetric model was developed to predict the critical conditions for spontaneous combustion in wood pellet storage structures. The sensitivity of self heating to variable boundary conditions and storage dimensions was analyzed. The following conclusions can be drawn from this study. 8.1.1. Measurement of thermal properties The results showed that moisture content has a significant effect on both effective thermal conductivity and specific heat capacity in the range of 1.4 to 9% w.b.. The effective thermal conductivity of wood pellets ranges from 0.139 to 0.192 W/(m∙K), and the specific heat capacity of wood pellets ranges from 1.074 to 1.253 kJ/(kg∙K), with both increasing with increasing pellet  moisture. The effective thermal conductivity is more sensitive to moisture than to porosity and insensitive to the size of pellets. Empirical correlations were developed based on experimental results to predict the effective thermal conductivity and specific heat capacity of pellets as a function of pellet moisture content and porosity. 8.1.2.  Exothermic properties of pellets Temperature and age of pellets are two major factors impacting self-heating phenomena. A higher heat release rate is associated with a higher internal temperature. Older pellets release less heat compared to freshly made pellets. A maximum heat release rate of 0.33 mW/g was observed at 50oC for fresh pellets. However, the heat release rate is less sensitive to pellet moisture. 159  The estimated activation energy of fresh wood pellets ranges from 50-60 kJ/mol. The activation energy also increases with the age of pellets. A kinetic equation for heat release rate was developed. This equation was implemented in simulation model to analyze self- heating at different storage temperatures and ages of wood pellets. 8.1.3. Wood pellets off-gassing during storage Experiments were carried out to study gas emissions during bulk storage of wood pellets at temperature 30oC to 50oC. The results show that the off-gassing from bulk storage of wood pellets mainly consists of CO2, CO, CH4 and H2. The gas emissions approximately follow a first order reaction kinetics with maximum concentrations of 2.3%, 1.88%, 0.193% and 1.5% (volumetric), respectively, and plateau emission factors of 0.0273, 0.0167, 0.0008 and 0.0007 g/kg (pellets), respectively. The emission factors of CO2, CO and CH4 increase with   reactor siz. Emission factor of H2 is less affected by the reactor size. The agre of Pellets   has significant impacts on CO2, CO, CH4 and H2 emissions. Emission factors of all gases decrease with increasing the age of pellets. However, the effect of pellets age decreases as the test temperature increases. 8.1.4. Thermal runaway and critical temperature Laboratory experiments were carried out to investigate the self-heating phenomenon at high temperatures (100oC to 200oC). When temperature reached 200oC, the heat generated from the pyrolysis reaction became significant enough to cause an invisible smoldering combustion within the bulk wood pellets. The minimum ambient (oven) temperature that could cause a thermal runaway inside of bulk wood pellets, defined as the critical ambient temperature, was found to be 135 oC, 160 oC, 173 oC and 191 oC for reactor sizes of radius 0.278, 0.108, 0.082 and 0.053 m, respectively. Thermal runaway temperature profiles were analyzed using, FK (Frank-Kamenetskii's) method and numerical curve-fitting method. Mean values for activation energy (E) of 78.7±0.8 kJ/mol and self-heating rate constant (∆rhA) of (4.22±2.5)*10 6 kJ/(kg s) were obtained for four white pellets samples from different pellet producers in British Columbia. 160  8.1.5. Simulation of self-heating in wood pellet silos Using experimental results, a two dimensional axi-symmetric model was developed to simulate the self-heating process and to predict the thermal runaway (spontaneous combustion) during wood pellet long-term storage. A parametric study was carried out to investigate the effects of airflow rate, wall insulation, wind speed, ambient temperature and storage silo dimension. The following conclusions are derived. o  Both wall thickness and convection on the outer wall have a minimal effect on thermal runaway. Wall thickness and convection outside of the storage wall only affect the temperature of locations near the wall, but does have no effect on inner locations, especially at the center. o  Dimensions of the storage container have a significant impact on thermal runaway and spontaneous combustion. The critical ambient temperatures for silos with a radius of 0.2-10 m showed that the critical ambient temperature for a 10 m radius silo could be as low as 36 oC. o  The inner air ventilation effect is investigated for both continuous airflow and periodical airflow. A critical air velocity is calculated for each given ambient temperature. The results showed that the air ventilation is a very effective technique for reducing self-heating and preventing thermal runaway at ambient temperatures lower than 330 K. The ventilation effect becomes very limited when the ambient temperature is higher than 330 K. o  The temperature profiles data from an industrial silo were collected to compare with the simulation results. The simulated temperature developing process from the model is very close to the industrial data. However, the airflow cooling process is faster and more effective in the simulation than observed in the industrial data.  However, due to the limitation of the experimental data, there are some limitations regarding the above conclusions. o The modeling conclusion may not apply at a higher humidity condition or when pellets are wet. Pellets trend to have a higher self-heating speed at higher humidity condition. 161  o The results are based on white pellets produced in British Columbia, Canada. The properties of pellets produced from other feedstock, such as bark may be different. o The model can accurately predict the self-heating rate and temperature development before the thermal runaway. However, after the thermal runaway the reaction rate and heating behavior may not be able to predict accurately due to the changes of reaction mechanism and the thermal properties of the materials. 162  8.2. Recommendations for industrial practitioners Based on the simulation results of this study, the following recommendations are made for the industrial practitioners. o The self-heating process needs to be simulated when designing a wood pellets storage system. Whether or not a thermal runaway will occur and the time to initiate a thermal runaway need to be determined for the system, and the ambient condition to prevent the dangerous phenomenon needs to be known. o Ventilation inside of the silo is essential when designing a wood pellets storage system. The cooling of the ventilation is significantly effective in lowering the temperature of the silo. However, the effective superficial air velocity should be controlled to be less than 0.5 m/s according to the simulation results of this study, because the temperature rise becomes insensitive to the ventilation rate beyond 0.5 m/s. o Increasing the pellets' age before storage will be one optional approach to reduce self-heating potential and to lower the pellet temperature at storage. o Since the self-heating is proportional to the silo size, it is recommended to use several smaller silos rather than one large silo in designing the storage system. 8.3. Recommendations for further work This research experimentally measured the self-heating during wood pellets storage and numerically simulated the self-heating developing process in wood pellets storage. During this study, effect of moisture content, bulk density, pellet age, pellet size are investigated. Self-heating of wood pellets is a  complex process and affected by multiple  factors. Further research  needs to be done to quantify all factors affecting wood pellets self heating . The following  is a set of  recommended for   future study. o  In the low temperature self-heating measurement, effect of relative humidity in the container has not been studied in this work. The relative humidity of air is an important factor affecting the self-heating, which needs to be investigated. 163  o  Because the maximum operation temperature of the TAM air calorimeter is 50oC, self-heating behavior at temperature range of 50oC to 100oC has not been investigated in this study, for example, the effect of moisture content and pellets age on the self-heating rate at this temperature range. More work need to be done. o  In the real commercial silo, the airflow inside the silo is not axi-symmetric as we assumed in the current simulation. A three-dimensional self-heating model needs to be developed in order to improve the accuracy of prediction of self-heating in commercial silos.  164    References 1. BC Ministry of Energy, Mines and Petroleum Resources, and BC Ministry of Forests and Range. An information guide on pursuing biomass energy opportunities and technologies in British Columbia. in Biomass & Bioenergy Conference2008. Prince Geroge. 2. Explosion at Armstrong plant, in www.globaltvbc.com. 2011: Kelowna. 3. Fire destroys wood stove business in Norwich 2011. [Acessed 27th Apr. 2012]. http://www.forestbusinessnetwork.com/10557/fire-destroys-wood-stove-business- in-norwich/. 4. Firefighters battle huge biomass fire at Port of Tyne 2011. [Acessed 27th Apr. 2012]. http://www.forestbusinessnetwork.com/10442/firefighters-battle-huge-biomass- fire-at-port-of-tyne/. 5. Explosion damages Waycross plant; no injuries reported. 2011. [Acessed 27th Apr. 2012]. http://www.forestbusinessnetwork.com/4314/explosion-damages-waycross-plant- no-injuries-reported/. 6. Firefighters leave power station, in BBC. 2012. 7. Power station engulfed in blaze, in BBC. 2012. 8. Firefighters battle blaze at Scotland plant 2012. [Acessed 27th Apr. 2012]. http://www.laurinburgexchange.com/view/full_story/17804602/article- Firefighters-battle-blaze-at-Scotland-plant?instance=popular. 9. Units fight another blaze at NC wood pellet plant, in WMBF News. 2012: North Carolina. 10. Wood Pellets Association of Canada, Material safety data sheet for wood pellets in bags. 2010, Shaw Resources: Canada. 11. Rowell, R.M., The chemistry of solid wood. Advances in chemistry series. 1983, Washington, D.C.: American Chemical Society. 12. Klemm, D., Heublein, B., Fink, H.-P., and Bohn, A., Cellulose: fascinating biopolymer and sustainable raw material. Angewandte Chemie International Edition, 2005. 44(22): p. 3358-3393. 13. Laine, C., Structures of hemicelluloses and pectins in wood and pulp. 2005, Helsinki Univeristy of Technology: Espoo, Finland. 14. Eriksson, Ö., Goring, D.A.I., and Lindgren, B.O., Structural studies on the chemical bonds between lignins and carbohydrates in spruce wood. Wood Science and Technology, 1980. 14(4): p. 267-279. 165  15. Teixeira, G., Van de Steene, L., Martin, E., Gelix, F., and Salvador, S., Gasification of char from wood pellets and from wood chips: Textural properties and thermochemical conversion along a continuous fixed bed. Fuel, 2012. 102(0): p. 514-524. 16. Lam, P.S., Steam explosion of biomass to produce durable wood pellets, Ph.D Thesis, University of British Columbia. 2011. 17. Lam, P.S., Sokhansanj, S., Bi, X., Lim, C.J., and Melin, S., Energy Input and Quality of Pellets Made from Steam-Exploded Douglas Fir (Pseudotsuga menziesii). Energy & Fuels, 2011. 25(4): p. 1521-1528. 18. Lavery, Measurement of VOC emissions from ponderosa pine lumber using commercial and laboratory kilns. Drying technology, 2001. 19(9): p. 2151. 19. Meijer, R. and Gast, C.H. Spontaneous combustion of biomass: experimental study into guidelines to avoid and control this phenomena. in 2nd World Conference on Biomass for Energy, Industry and Climate Protection.10-14 May, 2004. Rome, Italy 20. Arshadi, M. and Gref, R., Emission of volatile organic compounds from softwood pellets during storage. Forest Products Journal, 2005. 55(12): p. 132-135. 21. Wihersaari, M., Evaluation of greenhouse gas emission risks from storage of wood residue. Biomass and Bioenergy, 2005. 28(5): p. 444-453. 22. Larsson, S.H., Lestander, T.A., Crompton, D., Melin, S., and Sokhansanj, S., Temperature patterns in large scale wood pellet silo storage. Applied Energy, 2012. 92(0): p. 322-327. 23. Maheshwari, R., Bharadwaj, G., and Bhat, M.K., Thermophilic fungi: their physiology and enzymes. Microbiology and Molecular Biology Reviews, 2000. 64(3): p. 461-488. 24. Kubler, H., Heat generating processes as cause of spontaneous ignition in forest products. Forest Products Abstracts, 1987. 10: p. 299-327. 25. Rupar-Gadd, K., Biomass pre-treatment for the production of sustainable energy - emissions and self-ignition, Ph.D Thesis, Vaxjo University, Sweden. 2006. 26. Wadso, L., Measuring chemical heat production rates of biofuels by isothermal calorimetry for hazardous evaluation modelling. Fire and Materials, 2007. 31(4): p. 241-255. 27. Chico-Santamarta, L., Humphries, A.C., Chaney, K., White, D.R., Magan, N., and Godwin, R.J., Microbial changes during the on-farm storage of canola (oilseed rape) straw bales and pellets. Biomass and Bioenergy, 2011. 35(7): p. 2939-2949. 28. Lehtikangas, P., Storage effects on pelletised sawdust, logging residues and bark. Biomass and Bioenergy, 2000. 19(5): p. 287-293. 29. Yazdanpanah, F., Sokhansanj, S., Lau, A.K., Lim, C.J., Bi, X., and Melin, S., Airflow versus pressure drop for bulk wood pellets. Biomass and Bioenergy, 2011. 35(5): p. 1960- 1966. 166  30. Yazdanpanah, F., Sokhansanj, S., Lau, A.K., Lim, C.J., Bi, X., Melin, S., and Afzal, M., Permeability of wood pellets in the presence of fines. Bioresource Technology, 2010. 101(14): p. 5565-5570. 31. Pa, A. and Bi, X.T., Modeling of off-gas emissions from wood pellets during marine transportation. Annals of Occupational Hygiene, 2010. 54(7): p. 833-841. 32. Svedberg, U., Samuelsson, J., and Melin, S., Hazardous off-gassing of carbon monoxide and oxygen depletion during ocean transportation of wood pellets. Annals of Occupational Hygiene, 2008. 52(4): p. 259-266. 33. Kuang, X., Shankar, T.J., Bi, X.T., Lim, C.J., Sokhansanj, S., and Melin, S., Rate and peak concentrations of off-gas emissions in stored wood pellets—sensitivities to temperature, relative humidity, and headspace volume. Annals of Occupational Hygiene, 2009. 53(8): p. 789- 796. 34. Kuang, X., Shankar, T.J., Bi, X.T., Sokhansanj, S., Jim Lim, C., and Melin, S., Characterization and kinetics study of off-gas emissions from stored wood pellets. Annals of Occupational Hygiene, 2008. 52(8): p. 675-683. 35. Kuang, X., Shankar, T.J., Sokhansanj, S., Lim, C.J., Bi, X.T., and Melin, S., Effects of headspace and oxygen level on off-gas emissions from wood pellets in storage. Annals of Occupational Hygiene, 2009. 53(8): p. 807-813. 36. Svedberg, U.R.A., Hogberg, H.E., Hogberg, J., and Galle, B., Emission of hexanal and carbon monoxide from storage of wood pellets, a potential occupational and domestic health hazard. Annals of Occupational Hygiene, 2004. 48(4): p. 339-349. 37. Risholm-Sundman, M.L., M.; Vestin, E.; Herder, P., Emissions of acetic acid and other volatile organic compounds from different species of solid wood. Holz als Roh- und Werkstoff, 1998. 56: p. 125-129. 38. Wihersaari, M., Greenhouse gas emissions from final harvest fuel chip production in Finland. Biomass and Bioenergy, 2005. 28(5): p. 435-443. 39. Manninen, A.-M., Pasanen, P., and Holopainen, J.K., Comparing the VOC emissions between air-dried and heat-treated Scots pine wood. Atmospheric Environment, 2002. 36(11): p. 1763-1768. 40. Banerjee, S., Mechanisms of terpene release during sawdust and flake drying. Holzforschung, 2001. 55(4): p. 413-416. 41. Martínez-Inigo, M.J., Immerzeel, P., Gutierrez, A., Río, J.C.d., and Sierra-Alvarez, R., Biodegradability of Extractives in Sapwood and Heartwood from Scots Pine by Sapstain and White-Rot Fungi. Holzforschung, 1999. 53(3): p. 247-252. 42. Ryckeboer, J., Mergaert, J., Vaes, K., Klammer, S., De Clercq, D., Coosemans, J., Insam, H., and Swings, J., A survey of bacteria and fungi occurring during composting and self- heating processes. Annals of Microbiology, 2003. 53(4): p. 349-410. 167  43. Petrucci, R.H., Hardwood, W.S., Herring, F.G., and Madura, J.D., Gerneral chemistry. 9th ed. Principles and Modern Applications. 2007, New Jersey: Pearson Education, Inc. 44. Richardson, J.B., R; et al, Bioenergy from sustainable Forestry. 2002, Dordrecht/ Boston/ London: Kluwer Academic Publishers. 45. Yang, L.Z.G.Y.Z.W.F., The influence of different external heating ways on pyrolysis and spontaneous ignition of some woods. Journal of Analytical and Applied Pyrolysis, 2007. 78(1): p. 40-45. 46. Fasina, O.O. and Sokhansanj, S., Bulk thermal properties of alfalfa pellets. Canadian Agricultural Engineering, 1995. 37(2): p. 91. 47. Lynd, L.R., Elander, R.T., and Wyman, C.E., Likely features and costs of mature biomass ethanol technology. Applied Biochemistry and Biotechnology, 1996. 57-8: p. 741-761. 48. Lynd, L.R., Overview and evaluation of fuel ethanol from cellulosic biomass: Technology, economics, the environment, and policy. Annual Review of Energy and the Environment, 1996. 21: p. 403-465. 49. Kobayashi, T. and Sakai, Y., Hydrolysis rate of pentosan of hardwood in dilute sulfuric acid. Bulletin of the Agricultural Chemical Society of Japan, 1956. 20(1): p. 1-7. 50. Saeman, J.F., Kinetics of wood saccharifica Kinetics of wood saccharification - hydrolysis of cellulose and decomposition of sugars in dilute acid at high temperature. Industrial and Engineering Chemistry, 1945. 37(1): p. 43-52. 51. Burns, D.S., Ooshima, H., and Converse, A.O., Surface-area of pretreated lignocellulosics as a function of the extent. Applied Biochemistry and Biotechnology, 1989. 20-1: p. 79-94. 52. Young, R.A. and Rowell, R.M., Cellulose: structure, modification, and hydrolysis. 1945, New York: Wiley. 53. Jacobsen, S.E. and Wyman, C.E., Cellulose and hemicellulose hydrolysis models for application to current and novel pretreatment processes. Applied Biochemistry and Biotechnology, 2000. 84-6: p. 81-96. 54. Mohan, D., Pittman, C.U., and Steele, P.H., Pyrolysis of wood/biomass for bio-oil: A critical review. Energy & Fuels, 2006. 20(3): p. 848-889. 55. Yang, H., Yan, R., Chen, H., Lee, D.H., and Zheng, C., Characteristics of hemicellulose, cellulose and lignin pyrolysis. Fuel, 2007. 86(12-13): p. 1781-1788. 56. Bridgwater, A.V. and Peacocke, G.V.C., Fast pyrolysis processes for biomass. Renewable & Sustainable Energy Reviews, 2000. 4(1): p. 1-73. 57. Bilbao, A model for the prediction of the thermal degradation and ignition of wood under constant and variable heat flux. Journal of Analytical and Applied Pyrolysis, 2002. 62(1): p. 63. 168  58. Bilbao, R., Mastral, J.F., Aldea, M.E., Ceamanos, J., Betr 醤, M., and Lana, J.A., Experimental and theoretical study of the ignition and smoldering of wood including convective effects. Combustion and Flame, 2001. 126(1-2): p. 1363-1372. 59. Manya, J.J., Velo, E., and Puigjaner, L., Kinetics of biomass pyrolysis: A reformulated three- parallel-reactions model. Industrial & Engineering Chemistry Research, 2003. 42(3): p. 434-441. 60. Varhegyi, G., Antal, M.J., Jakab, E., and Szabo, P., Kinetic modeling of biomass pyrolysis. Journal of Analytical and Applied Pyrolysis, 1997. 42(1): p. 73-87. 61. Shafizadeh, F., The Chemistry of Pyrolysis and Combustion, in The Chemistry of Solid Wood. 1984, American Chemical Society. p. 489-529. 62. Ohlemiller, T.J. Smoldering combustion propagation on solid wood. in Fire Safety Science. Proceedings. 3rd International Symposium1991. Edinburgh, Scotland: International Association for Fire Safety Science. 63. Ohlemiller, T.J., Modeling of smoldering combustion propagation. Progress in Energy and Combustion Science, 1985. 11(4): p. 277-310. 64. Buckmaster, J. and Lozinski, D., An elementary discussion of forward smoldering. Combustion and Flame, 1996. 104(3): p. 300-310. 65. McCarter, R.J., Smoldering combustion of wood fibers - cause and prevention. Journal of Fire & Flammability, 1978. 9(1): p. 119-126. 66. Miao, Mechanism of spontaneous heating of hay .1. Necessary conditions and heat-generation from chemical-reactions. Transactions of the ASAE, 1994. 37(5): p. 1561. 67. Liodakis, The effect of (NH4)(2)HPO4 and (NH4)(2)SO4 on the spontaneous ignition properties of Pinus halepensis pine needles. Fire safety journal, 2002. 37(5): p. 481. 68. ASTM Standard, E1491-06 Standard test method for minimum autoignition temperature of dust clouds. 2006: ASTM International. 69. Melin, S., Testing of explosibility and fammability of airborne dust from wood pellets. 2008, Wood pellet association of canada. 70. Wood pellet association of canada, Material safety data sheet of wood pellets in bulk. 2011. 71. Calle, S., Klaba, L., Thomas, D., Perrin, L., and Dufaud, O., Influence of the size distribution and concentration on wood dust explosion: Experiments and reaction modelling. Powder Technology, 2005. 157(1-3): p. 144-148. 72. Soundararajan, R., Amyotte, P.R., and Pegg, M.J., Explosibility hazard of iron sulphide dusts as a function of particle size. Journal of Hazardous Materials, 1996. 51(1–3): p. 225- 239. 73. Cashdollar, K.L., Coal dust explosibility. Journal of Loss Prevention in the Process Industries, 1996. 9(1): p. 65-76. 169  74. S.A. Sawmilling Pty Ltd, Material safety data sheet of sawdust from non-treated softwood. . 2000. 75. Hyne & Son PTY Limited, Material safety data sheet of sawdust from non treated hardwood, in Fuel. 2007: Australia. 76. Yan, Z.H. and Nilsson, E.E.A., Large eddy simulation of natural convection along a vertical isothermal surface. Heat and Mass Transfer, 2005. 41(11): p. 1004-1013. 77. Nelson, M.I., Balakrishnan, E., and Chen, X.D., A Semenov Model of Self-Heating in Compost Piles. Process Safety and Environmental Protection, 2003. 81(5): p. 375-383. 78. Bilbao, R., Millera, A., and Arauzo, J., Kinetics of weight loss by thermal decomposition of different lignocellulosic materials. Relation between the results obtained from isothermal and dynamic experiments. Thermochimica Acta, 1990. 165(1): p. 103-112. 79. Bilbao, R., Millera, A., and Murillo, M.B., Temperature profiles and weight loss in the thermal decomposition of large spherical wood particles. Industrial & Engineering Chemistry Research, 1993. 32(9): p. 1811-1817. 80. Scott., F.H., Elements of chemical reaction engineering. 4th Edition ed. 2005, Ann Arbor: Pearson Education, Inc. 81. Tzeng, L.S., Atreya, A., and Wichman, I.S., A one-dimensional model of piloted ignition. Combustion and Flame, 1990. 80(1): p. 94-107. 82. Mohsenin, N.N., Thermal properties of foods and agricultural materials. 1980: Gordon and Breach. 83. Adl-Zarrabi, B., Boström, L., and Wickström, U., Using the TPS method for determining the thermal properties of concrete and wood at elevated temperature. Fire and Materials, 2006. 30(5): p. 359-369. 84. Blomquist, P., Spontaneous ignition of biofuels-an experimental investigation through small- and large-scale tests. 2006, SP Technical Research Institue of Sweden. 85. Kollman, F. and Cote, W., Principles of wood science and technology. 1968, Berlin, Germany: Springer-Verlag. 86. Gupta, M., Yang, J., and Roy, C., Specific heat and thermal conductivity of softwood bark and softwood char particles. Fuel, 2003. 82(8): p. 919-927. 87. Gupta, M., Yang, J., and Roy, C., Predicting the effective thermal conductivity of polydispersed beds of softwood bark and softwood char. Fuel, 2003. 82(4): p. 395-404. 88. Peters, B. and Bruch, C., Drying and pyrolysis of wood particles: experiments and simulation. Journal of Analytical and Applied Pyrolysis, 2003. 70(2): p. 233-250. 89. Vijeu, R., Gerun, L., Tazerout, M., Castelain, C., and Bellettre, J., Dimensional modelling of wood pyrolysis using a nodal approach. Fuel, 2008. 87(15-16): p. 3292-3303. 170  90. Pauner, M. and Bygbjerg, H., Spontaneous ignition in storage and production lines: Investigation on wood pellets and protein powders. Fire and Materials, 2007. 31(8): p. 477-494. 91. Blanchard, E., Simulation of self-heating in wood pellet storage using input data from small-scale experiments. 2007, SP Technical Research Institute of Sweden. 92. Walker, I. and Harrison, W., Self-heating of wet wood. 1. Exothermic oxidation of wet sawdust. New Zealand Journal of Science, 1977. 20: p. 191-200. 93. Bunyan, P., Baker, C., and Turner, N., Application of heat conduction calorimetry to high explosives. 3rd Internation Symposium on the Heat Flow Calorimetry of Energetic Materials, 2003. 401(1): p. 9-16. 94. Willson, R.J., Beezer, A.E., Mitchell, J.C., and Loh, W., Determination of thermodynamic and kinetic parameters from isothermal heat conduction microcalorimetry: applications to long-term- reaction studies. The Journal of Physical Chemistry, 1995. 99(18): p. 7108-7113. 95. Jones, J.C., Chiz, P.S., Koh, R., and Matthew, J., Kinetic parameters of oxidation of bituminous coals from heat-release rate measurements. Fuel, 1996. 75(15): p. 1755-1757. 96. Jones, J.C., Towards an alternative criterion for the shipping safety of activated carbons. Journal of Loss Prevention in the Process Industries, 1998. 11(6): p. 407-411. 97. Bunyan, P.F., Griffiths, T.T., and Norris, V.J., Combined use of adiabatic calorimetry and heat conduction calorimetry for quantifying propellant cook-off hazards. 3rd Internation Symposium on the Heat Flow Calorimetry of Energetic Materials, 2003. 401(1): p. 17-25. 98. Hofelich, T.C. and LaBarge, M.S., On the use and misuse of detected onset temperature of calorimetric experiments for reactive chemicals. Journal of Loss Prevention in the Process Industries, 2002. 15(3): p. 163-168. 99. TA Instruments, TAM air isothermal calorimetry manual. 2011. 100. Cole, A., Chip pile fire. Pulp Paper Magazine of Canada, 1972. 73(C): p. 137-140. 101. Bowes PC, Selfheating: evalution and cotrolling hazards. Vol. 2. 1984, Quincy, MA: Dept. of the Environment, Building Research Establishment ; Amsterdam ; New York : Elsevier. 180-189. 102. Frank-Kamenetskii, D., Diffusion and heat transfter in chemical kinetics (2nd edn) 1969, New York. 103. Beever PF, SFPE Handbook of fire protection engineering (2nd edn). Self heating and spontaneous combustion. 1988. 180-189. 104. Pauner, M.B.H., Spontaneous ignition in storge and production lines, Part 6: investagation of 6 mm and 8 mm wood pellets. DIFT Report 2006:02, Danish Institute of Fire and Security Technology, http://www.dift.dk/Research_reports_2006.asp, 2006. 171  105. Pauner, M.B.H., Spontaneous ignition in storge and production lines, Part 5: investagation of 8 mm wood pellets. DIFT Report 2006:01, Danish Institute of Fire and Security Technology, http://www.dift.dk/Research_reports_2006.asp, 2006. 106. Pauner, M.B.H., Spontaneous ignition in storge and production lines, Part 4: investagation of 6 mm wood pellets DIFT Report 2005:06, Danish Institute of Fire and Security Technology, http://www.dift.dk/Research_reports_2005.asp, 2005. 107. ASTM Standards, Standard test method for measurement of heat of hydration of hydraulic cementitious materials using isothermal conduction calorimetry. 2010. 108. TA Instrument, TAM Air Isothermal Calorimeter Operation Manual. 109. ASABE S269.4, Cubes, pellets and crumbles—definitions and methods for determining density, durability, and moisture content. 2007. 110. ASABE S358, S., Moisture Measurement—Foages. 2007. 111. ASABE Standards, Moisture measurement-foages. 2008. 112. Cengel, Y.A., Heat and mass transfer. 3rd ed. A practical approach. 2007: McGraw-Hill. 113. Tremblay, C., Cloutier, A., and Fortin, Y., Experimental determination of the convective heat and mass transfer coefficients for wood drying. Wood Science and Technology, 2000. 34(3): p. 253-276. 114. Bouguerra, A., Temperature and moisture dependence on the thermal conductivity of wood-cement- based composite: experimental and theoretical analysis. Journal of Physics D: Applied Physics, 1999. 32(21): p. 2797. 115. Cote, J. and Konrad, J.M., Assessment of structure effects on the thermal conductivity of two- phase porous geomaterials. International Journal of Heat and Mass Transfer, 2009. 52(3- 4): p. 796-804. 116. Kunii, D. and Smith, J.M., Heat transfer characteristics of porous rocks. AIChE Journal, 1960. 6(1): p. 71-78. 117. Woodside, W. and Messmer, J.H., Thermal conductivity of porous media. Part 1. unconsolidated sands. Journal of Applied Physics, 1961. 32(9): p. 1688-1699. 118. Deissler, R. and Boegli, J., An investigation of effective thermal conductivities of powders in various gases. ASME Transaction, 1958. 80: p. 1417-1425. 119. Tumuluru, J.S., Sokhansanj, S., Lim, C.J., Bi, T., Lau, A., Melin, S., Sowlati, T., and Oveisi, E., Quality of wood pellets produced in British Columbia for export. Applied Engineering in Agriculture, 2010. 26(6): p. 1013-1020. 120. Ingersoll, L.R., Zobel, O.J., and Ingersoll., A.C., Heat conduction with engineering and geological applications. 1950, New York: McGraw Hill Publisher. 121. Hooper, F.C. and Lepper, A.R., Transient heat flow apparatus for the determination of thermal conductivity. ASHVE Transactions, 1950. 56: p. 309-324. 172  122. Yang, W. and Sokhansanj, S., Determination of thermal conductivity, specific heat and thermal diffusivity of borage seeds. Biosystems Engineering, 2002. 82(2): p. 169. 123. Lagarias, J.C., Reeds, J.A., Wright, M.H., and Wright, P.E., Convergence properties of the Nelder--Mead simplex method in low dimensions. SIAM Journal on Optimization, 1998. 9(1): p. 112-147. 124. Bennett, C.A. and Franklin, N.L., Statistical analysis in chemistry and the chemical industry. 1954, New York: John Wiley & Sons. Inc. . 125. Guo, W., Lim, C.J., Bi, X., Sokhansanj, S., and Melin, S., Determination of effective thermal conductivity and specific heat of wood pellets. Fuel, 2013. 103: p. 347-355. 126. Jaya Shankar, T., Kuang, X., Shahab, S., Lim, C.J., Tony, B., and Staffan, M., Effect of storage temperature on off-gassing and physical properties of wood pellets, in 2008 Providence, Rhode Island, June 29 – July 2, 2008. 2008. 127. Jaya Shankar, T., Kuang, X., Shahab, S., Lim, C.J., Tony, B., and Staffan, M., Wood pellet properties produced in British Columbia and studies on off-gas emissions at low and high storage temperatures, in 2009 ASABE Aunal International Meeting, June 21 - June 24. 2009: Reno, Nevada. 128. Wark, K., Warner CF., Davis WT., Air pollution: its origin and control. 3rd edn. ed. 1999, Berkeley CA: Addison-Wesley. 337. 129. Guo, W., Trischuk, K., Bi, X., Lim, C.J., and Sokhansanj, S., Measurements of wood pellet self-heating kinetic parameters using isothermal calorimetry. Biomass Bionerergy, 2011. Under Review. 130. COMSOL 4.1, COMSOL Mutiphsics User Guide. 2011, COMSOL Inc. . 131. Bird, R.B., Transport phenomena 2nd Edition ed. 2002: Edwin N. Lightfoot. 282. 132. Yazdanpanah, F., Sokhansanj, S., Lim, C.J., Lau, A., Bi, X., and Melin, S., Stratification of off-gases in stored wood pellets. Biomass & Bioenergy, 2012. Under review. 133. Shankar, T.J., Sokhansanj, S., Lim, J., Bi, T., Melin, S., and Swaan, J., Design and instrumentation of a large column container to study the off-gassing, gas stratification and physical properties of wood pellets during storage, in ASABE Annual International Meeting. 2009: Reno, Nevada.   173  Appendix A Model parameters A.1. Material properties Wood pellets Properties  Expressions / Values Thermal conductivity, λ, W/(mK) 1 (0.219 0.01 ) (1 ) 0.027M        Specific heat capacity, CP, J/(kgK) P,1C = (1- M /100)(0.00546T -0.524) + 0.042M Density, kg/m3 650-700 Self-heating rate, Q=Q1+Q2, J/(kgK) With enough oxygen 11 11637.80 ( ) 1.15 15exp( 0.0256 ) 0.0026Q T e x T        9 2 ( ) 4.22 10 exp[ 78700 / (8.314 )]Q T T     Self-heating rate, Q=Q1+Q2, J/(kgK) With enough oxygen 1 1 1 11637.80 ( ) 1.15 15exp( 0.0256 ) exp[ 2.883exp( 0.014 6801.78 / ) 1 4 )] 0.0026 Q T e x T x T e t               9 2 ( ) 4.22 10 exp[ 78700 / (8.314 )]Q T T     Average particle diameter, DS, m 0.0061 Average particle length, LS,m 0.010 ∅ 0.42 a 477  174  Air Properties Expressions Specific heat capacity, CP, J/(kgK) 0.0769 1076.9PC T  1.1e3 Density, kg/m3 328.8 / (8.314 )P e T    1.22 Thermal conductivity, 𝛌, W/(mK) 0.8616log10( ( )) 3.714210 abs T    0.0251 Dynamic Viscosity, 𝛍, Pa·s 12 2 8 67.887 10 4.427 10 5.204 10T T              1.7e-5 Steel and Concrete Properties  Steel Concrete Thermal conductivity, 𝛌, W/(mK) 44.5 1.8 Specific heat capacity, CP, J/(kgK) 475 750 Density, kg/m3 7850 2300  175  Appendix B Industrial silo information The wood pellet silos are belong to the Fiberco Export Inc. and located at 1209 McKeen Avenue North Vancouver, B.C., Canada, V7P 3H9. The temperatures inside of the silo are monitored constantly using 12 cables suspended from the top at different locations with sensors at 10 levels for each cable. Figure B.1 shows a diagram of one silo and the locations of the cables. The bottom of the silo is 21.9 m and the total height of the silo is 23.2 m. However, 23.2 m of total height is composed of 17.1 m of cylindrical height and 6.1 m of circular cone roof. In the simulation, only the 17.1 m of cylindrical is consider ed.   Figure B.1 Diagram of the silo dimension and cable locations. (Reprinted from Larsson et. al. [22] with permission from Elsevier) 176   Appendix C Pilot scale experiments C.1 Equipments A pilot experiment has been conducted in a cylindrical reactor of 1.2 m in diameter and 4.6 m in height in order to verify the simulation results. The reactor is equipped with 30 thermocouples located at 6 different levels as shown in Figure C.1. The thermocouples are K-type thermocouples and their tips located at T0-T29 as indicated in Figure C.1. The temperatures of T0-T29 were recorded using a self-coded computer program with frequency of 1 minute. At the beginning of the experiment, 3 metric tons fresh wood pellets were loaded from the top. The reactor was sealed after loading and temperatures were monitored for more than 60 days. The experimental results of all thermocouples are shown in Figure C.2 and Figure C.3. During the experiments, the off-gassing inside of the silo was also measure continuously by my colleague, Fahimeh Yazdanpanah[132]. The O2 concentration of two levels is shown in Figure C.4. 177   Figure C.1 Diagram of the pilot size reactor and the locations of the thermocouples[133]. 178  0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level0left   Center T2 Half T3 Wall T4 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level1left   Center T7 Half T8 Wall T9 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level2left   Center T12 Half T13 Wall T14 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level3left   Center T17 Half T18 Wall T19 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level4left   Center T22 Half T23 Wall T24 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level5left   Center T27 Half T28 Wall T29  Figure C.2 Temperature profile from experimental results of lab scale experiments (Left side thermocouples) 179  0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level0right   Center T2 Half T1 Wall T0 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level1right   Center T7 Half T6 Wall T5 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level2right   Center T12 Half T11 Wall T10 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level3right   Center T17 Half T16 Wall T15 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level4right   Center T22 Half T21 Wall T20 0 100 200 300 400 500 18 20 22 24 26 28 30 Time (hour) T e m p e ra tu re  ( o C ) level5right   Center T27 Half T26 Wall T25  Figure C.3 Temperature profile from experimental results of lab scale experiments (Right side thermocouples) 180    Figure C.4 O2 concentration inside of the silo during the experiments. (Reprinted from Yazdanpanah et. al. [132] with permission from Yazdanpanah F.)  181   Appendix D Program This simulation was worked in COMSOL with MATLAB environment. By invoking commands and modules from COMSOL Multiphysics in the MATLAB interface,  the two dimensional axis-symmetric model is developed.  Here is an example script for simulating the temperature development of wood pellets self-heating inside of a 20 meter diameter silo with initial temperature, 'T_init' and inside airflow 'v1'. function femsol_airflow0902(T_init,cas1,v1) % COMSOL Multiphysics Model M-file flclear fem  % COMSOL version clear vrsn vrsn.name = 'COMSOL 3.4'; vrsn.ext = ''; vrsn.major = 0; vrsn.build = 248; vrsn.rcs = '$Name:  $'; vrsn.date = '$Date: 2007/10/10 16:07:51 $'; fem.version = vrsn;  saveid=0; T_inf=T_init; %---------------------save folder------------------------- cas=char(65+cas1); tinterv=3600; %tag=['Tinf' num2str(T_inf)]; ss='D:\SH Model\MATLAB'; if saveid==1 xlswrite([ss '\results' '\' datestr(date,'mmmdd') '.xls'], tinterv, tag, [cas num2str(1)]); end directory_name='E:\Comsol results'; resdir=[directory_name '\' datestr(date,'mmmdd')];  182  if exist(resdir,'dir') else     mkdir(directory_name,datestr(date,'mmmdd')); end  if exist(datestr(date,'mmmdd'),'dir') else     mkdir(datestr(date,'mmmdd')); end  %-----------------------hloc-------------------------------------- %calculate the heat transfer coefficient between wood pellets and air % v1 is the airflow velocity of air inside of silo hloc=37*v1^0.33+88*v1^0.619; hloc1=12.21*v1^(-0.67)+54.5*v1^(-0.318);  %-----------------------cal hwall------------------------- %calculate the heat transfer coefficient between outside wall and air ro=10;  % radius of the silo %T_init=300; v0=0.1; % wind speed outside of silo H1=2*ro; P=10^5; T1=T_inf; CP=0.0769*T1+1076.9; %[J/kg K] rho=P*28.8e-3/8.314/T1;  %[kg/m3] k=10^(0.8616*log10(abs(T1))-3.7142); %[W/(mK)] eta= (-7.887E-12*T1^2+4.427E-08*T1+5.204E-06); %[Pa•s] Gr=rho^2*9.18/T1*H1^3/eta^2;   %%%% need to time deltaT Re=H1*v0*rho/eta; Pr=CP*eta/k; Nu=0.3975*(Gr*Pr)^0.25; % Nu need to multiply deltaT^0.25 Nu2=0.6*Re^0.5*Pr^(1/3); h2=k/H1;  183  %------------------------ Geometry---------------------------- %define geometry g3=rect2('10','13','base','corner','pos',{'0','0'},'rot','0'); g5=rect2('0.2','13','base','corner','pos',{'10','0'},'rot','0');  % Analyzed geometry clear s s.objs={g3,g5}; s.name={'R1','R2'}; s.tags={'g3','g5'};  fem.draw=struct('s',s); fem.geom=geomcsg(fem);  % Initialize mesh fem.mesh=meshinit(fem, ...                   'hauto',5); %figure, meshplot(fem) % (Default values are not included)  % Application mode 1 clear appl appl.mode.class = 'HeatTransfer'; appl.mode.type = 'axi'; appl.assignsuffix = '_ht'; clear bnd bnd.type = {'q0','cont','ax','q'}; bnd.Tinf = {0,0,0,T_init}; bnd.h = {0,0,0,'hcov'}; %--------------------------------------------- bnd.ind = [3,1,1,2,1,1,4]; appl.bnd = bnd; clear equ equ.C = {'mat1_cp_fun(T)','mat2_C'}; equ.init = T_init; equ.k = {'mat1_k','mat2_k'}; 184  equ.Q = {'Q1+Q2-Q_cc',0}; equ.rho = {'mat1_rho','mat2_rho'}; equ.ind = [1,2]; appl.equ = equ; fem.appl[130] = appl;  % Application mode 2 clear appl appl.mode.class = 'ConvCond'; appl.mode.type = 'axi'; appl.dim = {'T2'}; appl.module = 'CHEM'; appl.assignsuffix = '_chcc'; clear prop clear weakconstr weakconstr.value = 'off'; weakconstr.dim = {'lm2'}; prop.weakconstr = weakconstr; appl.prop = prop; clear bnd bnd.type = {'q0','cont','ax','qc','T'}; bnd.T0 = {0,0,0,0,T_init}; bnd.ind = [3,5,4,1,2,2,2]; appl.bnd = bnd; clear equ equ.eta = {'mat3_eta(T2[1/K])[Pa*s]',0}; equ.C = {'mat3_Cp(T2[1/K])[J/(kg*K)]',1006}; equ.gamma = {'mat3_gamma',1}; equ.Q = {'Q_cc',0}; equ.rho = {'mat3_rho(p[1/Pa],T2[1/K])[kg/m^3]',1.205}; equ.fluidtype = {'gas','userdefined'}; equ.sdiff = 'off'; equ.init = {T_init,0}; equ.k = {'mat3_k(T2[1/K])[W/(m*K)]',0.025}; equ.v = {'v',0}; 185  equ.u = {'u',0}; equ.usage = {1,0}; equ.ind = [1,2]; appl.equ = equ; fem.appl[109] = appl;  % Application mode 3 clear appl appl.mode.class = 'NavierStokes'; appl.mode.type = 'axi'; appl.module = 'CHEM'; appl.gporder = {4,2}; appl.cporder = {2,1}; appl.assignsuffix = '_chns'; clear prop prop.brinkmandef='On'; clear weakconstr weakconstr.value = 'off'; weakconstr.dim = {'lm3','lm4','lm5','lm6','lm7','lm8','lm9','lm10'}; prop.weakconstr = weakconstr; appl.prop = prop; clear bnd bnd.symtype = {'sym','sym','ax','sym','sym'}; bnd.v0 = {0,0,0,'v1',0}; bnd.p0 = {0,0,0,0,'10^5'}; bnd.type = {'walltype','int','sym','inlet','outlet'}; bnd.walltype = {'slip','noslip','noslip','noslip','noslip'}; bnd.U0in = {1,1,1,'v1',1}; bnd.ind = [3,4,5,1,2,2,2]; appl.bnd = bnd; clear equ equ.eta = {'mat3_eta(T2[1/K])[Pa*s]',1}; equ.gporder = {{1;1;2}}; equ.epsilonp = {0.45,'-39.658*v1+23.773'}; equ.rho = {'mat3_rho(p[1/Pa],T2[1/K])[kg/m^3]',1}; 186  equ.cporder = {{1;1;2}}; equ.init = {{0;0;0;'10^5';-8;-8;0;0;0;0},{0;0;0;0; ...   -8;-8;0;0;0;0}}; equ.gamma = {'mat3_gamma',1}; equ.sigma = {'mat3_sigma',0}; equ.usage = {1,0}; equ.ind = [1,2]; appl.equ = equ; fem.appl{3} = appl; fem.sdim = {'r','z'}; fem.frame = {'ref'}; fem.border = 1; clear units; units.basesystem = 'SI'; fem.units = units; fem.const = {'xx0',x0,...         'Nu11',Nu,...     'Nu22', Nu2,...     'hco',h2 ,...     'T_init',T_init,...     'hl',hloc...     'hl1',hloc1,...     'v1',v1}; % Global expressions fem.globalexpr = {     'x2','0.00021*(T-278)^3*(t/3600/24)',...%pellets age     'xx','x2+xx0',...     'Q1','650*(1.15e15*exp(-0.0256*(xx)-11637/T)+0.0026)', ...% heat gerneration from pellets      'Q2','650*4.22e9*exp(-78700/8.314/T)', ...%heat gerneration from pellets     'hloc','(hl+hl1*(v-v1))*477',...% calculate the h of pellets and air   'Q_cc','hloc*(T-T2)'}; % heat exchange between pellets and air fem.expr = {    'hcov','Nu11*hco*(T-T_init)+Nu22*hco' ...    }; % h of wall and ambient air 187   % Library materials clear lib lib.mat[130].name='Wood Pellets'; lib.mat[130].varname='mat1'; lib.mat[130].variables.C='cp_fun(T)'; lib.mat[130].variables.rho='650'; lib.mat[130].variables.k='0.165'; clear fcns fcns[130].type='inline'; fcns[109].type='inline'; fcns[109].name='cp_fun(T)'; fcns[109].expr='5.46*T-524'; fcns[109].dexpr={'diff(5.46*T-524,T)'}; lib.mat[130].functions = fcns; lib.mat[109].name='Steel AISI 4340'; lib.mat[109].varname='mat2'; lib.mat[109].variables.nu='0.28'; lib.mat[109].variables.E='205e9[Pa]'; lib.mat[109].variables.mur='1'; lib.mat[109].variables.sigma='4.032e6[S/m]'; lib.mat[109].variables.epsilonr='1'; lib.mat[109].variables.alpha='12.3e-6[1/K]'; lib.mat[109].variables.C='475[J/(kg*K)]'; lib.mat[109].variables.rho='7850[kg/m^3]'; lib.mat[109].variables.k='44.5[W/(m*K)]'; lib.mat{3}.name='Air, 1 atm'; lib.mat{3}.varname='mat3'; lib.mat{3}.variables.nu0='nu0(T[1/K])[m^2/s]'; lib.mat{3}.variables.eta='eta(T[1/K])[Pa*s]'; lib.mat{3}.variables.gamma='1.4'; lib.mat{3}.variables.sigma='0[S/m]'; lib.mat{3}.variables.C='Cp(T[1/K])[J/(kg*K)]'; lib.mat{3}.variables.rho='rho(p[1/Pa],T[1/K])[kg/m^3]'; lib.mat{3}.variables.k='k(T[1/K])[W/(m*K)]'; 188  lib.mat{3}.variables.cs='cs(T[1/K])[m/s]'; clear fcns fcns[130].type='inline'; fcns[130].name='cs(T)'; fcns[130].expr='sqrt(1.4*287*T)'; fcns[130].dexpr={'diff(sqrt(1.4*287*T),T)'}; fcns[109].type='inline'; fcns[109].name='nu0(T)'; fcns[109].expr='(-7.887E-12*T^2+4.427E-08*T+5.204E-06)/(1.013e5*28.8e-3/8.314/T)'; fcns[109].dexpr={'diff((-7.887E-12*T^2+4.427E-08*T+5.204E-06)/(1.013e5*28.8e- 3/8.314/T),T)'}; fcns{3}.type='inline'; fcns{3}.name='Cp(T)'; fcns{3}.expr='0.0769*T+1076.9'; fcns{3}.dexpr={'diff(0.0769*T+1076.9,T)'}; fcns{4}.type='inline'; fcns{4}.name='rho(p,T)'; fcns{4}.expr='p*28.8e-3/8.314/T'; fcns{4}.dexpr={'diff(p*28.8e-3/8.314/T,p)','diff(p*28.8e-3/8.314/T,T)'}; fcns[110].type='inline'; fcns[110].name='eta(T)'; fcns[110].expr='-7.887E-12*T^2+4.427E-08*T+5.204E-06'; fcns[110].dexpr={'diff(-7.887E-12*T^2+4.427E-08*T+5.204E-06,T)'}; fcns[109].type='inline'; fcns[109].name='k(T)'; fcns[109].expr='10^(0.8616*log10(abs(T))-3.7142)'; fcns[109].dexpr={'diff(10^(0.8616*log10(abs(T))-3.7142),T)'}; lib.mat{3}.functions = fcns;  fem.lib = lib;  % ODE Settings clear ode clear units; units.basesystem = 'SI'; 189  ode.units = units; fem.ode=ode; % Multiphysics fem=multiphysics(fem);  % Extend mesh fem.xmesh=meshextend(fem); t_list1=0:3600:3600*24*100; % Solve problem'init',T_init,... i=1; fem.sol=femtime(fem, ...                 'solcomp',{'v','T','u','p','T2'}, ...                 'outcomp',{'v','T','u','p','T2'}, ...                 'tlist',t_list1, ...                 'tout','tlist'); Tmax(i)=postmax(fem,'T'); %------------------save-------------------------------------------  tnf=[resdir '\results' num2str(T_init) 'v_' num2str(v1) '.mat']; save(tnf,'-struct','fem');  % Save current fem structure for restart purposes fem0=fem; sizeid=size(fem.sol.tlist,2); fem=multiphysics(fem);  end     190  Appendix E. Self-heating with free convection During the self-heating simulation,  it is assumed that the free convection effect on self- heating is negligible and can be safely omitted. In this section, a case is calculated, which simulated the self-heating with the existance of the free convection of the air inside of wood pellet silo, to examine the effect of free convection on the temperature profiles without forced air flow. E.1 Boundary condition A two dimensional axi-symmetric computational domain is shown in Figure E.1. The cylinderical storage domain has a diameter (D=10 meter) and height (H=8 meter). The wall thickness is set as Lw=0.1 meter. Case C : Free Convection B. 1 B. 3 B. 2 B. 4 Boundary Condition for Solid: B.1 Axial symmetry at r=0. B.2 Newton cooling law Eq.(7.1) at boundary B.2, B.3, B.4 (outside walls) Boundary Condition for Air: B.1 Axial symmetry at r=0. B.2 No slip wall for boundary B.2, B.3 and B.4. R=D/2 L W H LW Figure E. 1 Diagram of computational domain and boundary condition E.2 Governing equation For the solid phase (Wood Pellets) 191   2 1 1 1 1 ,1 1 12 1 ( ( ) ) ( )P G H T T T C r Q T Q t r r r z               (7.3) where 1 is the bulk density of the pellets, kg/m3; ,1P C  is the specific heat capacity of the bulk wood pellets, J/(kg K); T1 is temperature of the wood pellets, K; 1  is the effective thermal conductivity of wood pellets, W/(m K); QG(T1) is the heat generation from the bulk pellets, W/m3; QH is the convective heat flux from wood pellets to airflow, W/m 3; r and z are the cylindrical coordinates, m; t is the time, s. For gas phase (Air) Energy 2 2 2 2 2 2 ,2 2 2 1 ( ) ( ( ) )P z H T T T T C r Q t z r r r z                  (7.4) Momentum 2 2 2 22 1 ( ) ( ( ) ) ( )z z z zz r g T T t z r r r z                            (E.1) Continuity 2 2 ( ) 0z t z          (7.6) where 2 is the bulk density of the air, kg/m3; ,2P C  is the specific heat capacity of air, J/(kg K); T2 is temperature of air, K; 2  is thermal conductivity of air, W/(m K); z is advection velocity of the gas phase in z direction, m/s; µ2 is the viscosity of air, (Pa s); QH is the convective heat flux from wood pellets to airflow, W/m3; r and z are the cylindrical coordinates, m; t is the time, s. ϕ is the porosity of the pellets, in fraction.   is the average thermal Expansion coefficient, 1/K. E.3 Governing equation COMSOL is used to slove model euqations with the method described in section 7.3.2. E.3 Results 192  Ambient temperature of 320 K and ambient wind speed of 1m/s is used in this case in order to compare to the similar case  shown in Figure 7.10 without free convection. Figure E.2 shows the solid phase temperature profiles at different stages. Comparing to Figure 7.10, the solid phase temperature profiles are the same for  the first and second stages (3 days and 10 days). This is because the temperature profiles inside of wood pellets silo developed very slowly and evenly at the first and second stage, which give a small bouyant force 2 ( )g T T   . The air velocity caused by the free convection is around 10 -6 m/s as shown in Figure E.3. Therefore, free convection doesn’t have significant effect on the solid phase temperature profiles at these two initial stages. However, at the third and forth stages, after the thermal runaway , the temperature profiles are  different from those shown in Figure 7.10. The temperature gradients are smaller than those shown in Figure 7.10. This is because initial large temperature gradients (Ce and De) after thermal runaway cause a larger bouyant force 2 ( )g T T    which promotes free convection, improves heat transfer and reduces the temperature different within the silo. r,meter z ,m e te r A e  Time=3days   0 5 0 2 4 6 300 350 400 450 500 r,meter z ,m e te r B e  Time=10days   0 5 0 2 4 6 300 350 400 450 500 r,meter z ,m e te r C e  Time=20days   0 5 0 2 4 6 300 350 400 450 500 r,meter z ,m e te r D e  Time=23.2days   0 5 0 2 4 6 300 350 400 450 500  Figure E.2 Solid phase temperature profiles at different stages. 193    Figure E.3 Air velocity caused by the free convection at the 11th day of self-heating.    

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0073583/manifest

Comment

Related Items