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.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

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  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 selfheating 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.  ii  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 offgassing 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  iii  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.  iv  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  v  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  vi  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  vii  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 viii  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  ix  Appendix C Pilot scale experiments.............................................................................................. 176 Appendix D Program ...................................................................................................................... 181 Appendix E. Self-heaitng with free convection ........................................................................... 190  x  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  xi  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  xii  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  xiii  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 23oC 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  xiv  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 22oC 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 22oC 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 22oC 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 22oC 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 xv  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 ln( C )  2ln(T0 / r ) vs 1/ T0 for pellets A and Pellets E using FK method. ....................................................................................................................................... 112 Figure 6.5 Linear regression of ln(dT / dt ) P vs 1/ TP 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)*106 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  xvi  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  xvii  Nomenclature A1  Heat releasing rate coefficient, J/s  AS  Particles outer surface area, m2/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/(m2 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/m3  Nu  Nusselt number  P  Pressure, Pa  Pr  Prandtl number.  Q0  Heat rate per unit length of heater, W/m  xviii  QWall  Heat flux from the silo wall to the environment, W/m3  QG  Self-heating rate from wood pellets, W/m3  QH  The convective heat from wood pellets to airflow, W/m3  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  xix  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  xx  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.  xxi  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 1  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. Selfheating 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 selfheating 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.  2  Table 1.1 List of wood pellets fires accidents over last year  Location  Details  Category  Date  Armstrong Pellets plant[2]  Blast followed by visible fire spreading through the storage  Possible dust explosion  Apr.2011  British Columbia  building.  Shur Fire Energy  wood pellets plants storage  Possible spontaneous  Oct, 2011  Norwich[3]  heavy smoke, no visible fire and charring wall  combustion  Port of Tyne [4]  Fire happened inside of the wood pellets silo and very deep seated.  Spontaneous combustion  Nov 2011  Georgia Biomass plant [5]  Described as 'flash-type explosion’, which happened during  Dust explosion  Jun 2011  Germany  operation  Tilbury power station[6, 7]  Huge blaze broke out at power station in a storage area containing  Dust explosion  Feb 2012  United Kingdom  4000 tonnes of wood pellets  Laurinburg Nature' Earth plant  Fire in wood pellet storage silo and second fire one week later in  Possible spontaneous  Mar 2012  North Carolina[8, 9]  the same silo.  combustion  South Shields, United Kingdom  3  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 carbon4  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  5  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  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  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  Ultimate analysis (wt% dry basis)  LHV (MJ kg-1 dry basis)  6  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  7  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 selfheating, 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 2liter 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,  8  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].  9  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.  10  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].  11  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]. 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)  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 13  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  14  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].  15  Table 1.4 Results from explosibility testing dust from white pellets and bark pellets. (Reprinted from Melin [69] with permission from the author)  16  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 17  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.  18  Table 1.5 Published data and relationships on specific heat and thermal conductivity Thermal Conductivity, Researchers  Method and Material  W/(m K) and temperature range  Specific Heat, kJ/(kg K) and temperature range  dry  0.0002 dry  0.024 Kollman et al. [85]  General wood  [ room temperature]     (1  0.00015( M1  M 2 )) 2  Adl-Zarrabi et al. [83]  TPS method, A piece of wood  1  CP ,dry  0.0046T  0.116 [0 oC -100 oC ]  [5% <M< 35% ]  0.55 W/(m K) for parallel to grain direction 0.11 W/(m K) for perpendicular to grain direction  1.07 kJ/(kg K) for dry wood 1.38 kJ/(kg K) for M.C. 9.5% Room temperature  [20, 110 and 150 oC] Modified Fitch Gupta et al. [86]  apparatus and DSC, Dried softwood particles  0.0986 W/(m K) [37  oC]  CP  0.00546T  0.524 [40-140oC]  L=1mm D=6.5 mm Fasina et al. [46] Peters et al. [88]  Alfalfa pellets M.C. 7.5 to 18.0%  eff  0.049  0.0082M [20 oC]  CP  1.083  0.089M  0.0021M 2 [20 oC]  0.2 W/(m K) for No reference or  Fir wood particles (4mm)  CP=1.733 kJ/(kg K)  method reported  0.35 W/(m K) for  CP=1.112 kJ/(kg K)  Spruce wood particles Pauner et al.  Reference to  [90]  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.  19  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  20  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 10 4 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  T  2T CP   2  QG 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).  21  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:  QG (T )   r h  A exp(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   2T  2    r h  A exp( E / RT ) x  (1.3)  By introducing a dimensionless parameter δ, Eq. (1.3) is converted to a dimensionless form,   2   e 2   (1.4)  in which the dimensionless parameter δ is defined as  E r hA 2  RT0  re RT02  E  (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). 22  2.5 2.0 1.5    C  1.0 0.5 0.0 0.0  0.5  1.0  1.5  2.0   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 ln( C )  2ln(T0 / r ) as y coordinate and 1/ T0 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,  2.228  Side 2a=2√6 r  23  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   CP  T t      r h  A exp( E / RTP )  (1.6)  P  The identical temperature between this two location points is defined as the crossing point (P). After TP and T / t P (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(dT / dt ) P versus ( 1/ TP ). 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]  24  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.  25  . 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 selfheating 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  26  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 selfheating 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.  27  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  28  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.  Fibre glass insulation  66mm  T1 69 mm  T2  T3  113mm  610 mm T0  E, Electric heater  Fibre glass insulation Data logging  Electric heater controller  Computer  310 mm  Figure 2.1 Experimental setup for the measurement of thermal properties of wood pellets  29  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  30  Rubber Gasket  610 mm  Sampling Port  285 mm 275 mm  310 mm  160 mm  Sampling Port 1  45 L Steel cylinder  15 L Plastic bucket  Sampling Port 2 135 mm  255 mm  a. 45 L  b. 15 L Sampling Port  Sampling Port  105 mm  10 L Plastic bucket  265 mm  270 mm  240 mm  2 L glass cylinder  203 mm  c. 10 L  d. 2 L  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  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  32  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.  33  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  34  temperature was measured following the trial-and-error method. The reactor dimensions and thermocouple locations are illustrated in Figure 2.6.  Thermocouple 1 Thermocouple 2  Data Logging interface  Wood pellets container  oven  Figure 2.5 Experimental setup for measuring the critical oven (environment) temperature of wood pellet.  35  T3  T4  T2  T1  417 mm  88 mm  200 mm  58mm  116mm  356 mm  Self-heating reactor  T1  T2 T1  T2  164mm  T1 57 mm  40 mm  114 mm  94 mm  101 mm  188 mm  224 mm  5.6 mm  T2 28 mm  105 mm  216 mm  Large reactor  Medium reactor  Small reactor  Figure 2.6 Reactor dimensions and thermocouple locations.  36  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.  37  a. THELCO Precision laboratory oven  b. LHU-113 temperature and humidity cabinet  (Thermo Fisher Scientific Inc., 308 Ridgefield  (ESPEC North America, Inc., 4141 Central Parkway  Court Asheville, NC 28806)  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  38  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  39  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 LHU113 (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  40  Size Distribution for Sample 1 (Normal)  Frequency (%)  20  15  10  5  0  2  4  6  8  10  12  14  16  18  20  22  24  26  28  30  Pellet Length (mm)  a. Pellets A  Size Distribution for Sample 2A (Moisture content 4.6%) 25  Frequency (%)  20  15  10  5  0 0  5  10  15  20  25  30  35  40  Pellet Length (mm)  b. Original pellets B1 (Moisture content 4.6%)  Size Distribution for Sample 2B (Moisture content 8.8%) 25  Frequency (%)  20  15  10  5  0 0  5  10  15  20  25  30  35  40  Pellet Length (mm)  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  41  Samples  Pellet A  Pellet B Pellet B1  Pellet B2  640-670  690  590-620  1159  1186  1168  Average diameter, mm  6  6  6  Average length, mm  10  13.1  13.6  4.6%  4.6%  8.8%  1.17  1.22  1.23  0.023-0.090  0.024 -0.094  0.024-0.094  Density, kg/m3 Specific density,  kg/m3  Moisture content (as received), wt %, wet based. Lc, mm Biot Number  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 42  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.  80  Temperature(oC)  70  T0  60  50  T1 T2  40  T3  30  20  0  0.5  1  1.5  2  2.5  Time(s)  3  3.5  4  4.5 4  x 10  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  43  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,  Bi   hLC  s  (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,  LC   V  ( D / 2)2 L DL   2 AS  DL  2 ( D / 2) 4L  2D  (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  T  2T 1 T bCP  eff ( 2  )  evap t r r r  (3.3)  44  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/m3) in the air can be calculated from the water saturated vapor pressure using the idea gas law.  mv   n Psat ( RH )  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/m3) over a time range of Δt can be calculated by multiplying the moisture loss and the vaporization enthalpy hfg (kJ/kg).  Evap  evap  t  mv  Wwt  h fg  (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 E vap 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 55 oC 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.  45  Table 3.2. The ratio of evaporation heat (Evap) to total energy change (ρCpΔT) Tair,  oC  Psat  *,  RH  Pa  mv,  hfg* , kJ/kg  Evap, kJ/m3  ρCpΔT  Eevp/ρCpΔT  mol/m3  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 P sat 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,  bCP  (1   ) sCP ,s   f CP , f  (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/m3) at 20 oC is much smaller than the density of solid (ρs=1159 kg/m3), 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].  46  for parallel distribution and for serial distribution  eff   f  (1   )s 1  eff    (3.7)   1  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,  CP  (1- M / 100)CP , dry  ( M / 100)CP , M  (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.  bCP  T t  T 2   eff (  r  2    1 T r r  )  (3.10)  47  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/m3; 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  T    exp( x2 ) dx  T0 2eff  x Q0  where Q0  2 r0q0 is the heating rate per unit length of the heater , w/m;    (3.11)  r 2 eff t /  CP  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  ( t)  1/2  [82, 120-122].  48  T  T0  Q0B 2  eff  Q0 2  eff  1 r 2  Q0  1/2  4  ln(t )  (3.12)  eff  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  r  r0 , q0   eff T / r  B.C.2  r  R , T / r  0  I. C.  t  0 , T  T0  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 secondorder 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].  Residue   nt 150   (T nt 1  C   Ta )2  (3.13)  49  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 Mean and Variance  Trial 1  Trial 2  Trial 3  Trial 4  λ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  (at 95% confidence limit)1  Note: 1 Mean and Variance were calculated base on equation x=average(x)+st n-1,α/sqrt(n) [124]  50  38 Experimental Results Numerical Results Analytical Results  36  T1  Temperature(oC)  34 32  T2 30 28 26  T3  24 22  0  1000  2000  3000  4000  5000  6000  7000  8000  9000  Time(s) 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.  51  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.  52  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)  λeff, W/(m K)  λeff, W/(m K)  λS, W/(m K)  Experiments  Kollman2  Fasina3  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    eff  dry  0.0002 dry  0.024  and    (1  0.00015M 2 ) . dry  3  Fasina [46] effective thermal conductivity is calculated using equation eff  4  λS parallel is the single pellet thermal conductivity calculated using Eq.(3.7).   0.049  0.0082M .  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)  λeff,W/(m K)  λeff,W/(m K)  λS, W/(m K)  Experiments  Kollman2  Fasina3  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    eff  dry  0.0002 dry  0.024  and    (1  0.00015M 2 ) . dry  3  Fasina [46] effective thermal conductivity is calculated using equation eff  4  λS parallel is the single pellet thermal conductivity calculated using Eq. (3.7).   0.049  0.0082M .  53  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.  s  0.219  0.01M  (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.  eff  (0.219  0.01M )  (1   )    0.027  (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.  54  Single pellet thermal conductivity (W/(mk))  Pellet A Pellet B Eq. (3.14) Parallel distribution model Serial distribution model  0.40  0.35  0.30  0.25  0.20  0.15 0  2  4  6  8  10  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 s  water M / 100  (1  M / 100)s , dry Serial distribution model was calculated using equation  1  s    M 100  water    1  M 100  s , dry  where water  0.626 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.  55  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 (37oC) 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  CP  (1  M / 100)CP , dry  0.042M .  (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 CP  1.01  0.032M .  (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  56  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 curvefitting 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.  57  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.  58  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  1 year  1 year  Oven  2  2  Freeze  dried  months  months  dried  dried  M.C. stands for moisture contents of the pellets, wet based.  a b  Oven  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 0 oC). 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 23 oC) 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  M.C. stands for moisture contents of the pellets, wet based.  a  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  59  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  q  qox  qnox  (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.  60  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  rA    dCO 2  kCO 2 dt  (4.2)  where rA is the reaction rate base on O2 at time t, mol/( m3.s); k is the reaction rate constant, in s-1; CO2 is the concentration of O2 at time t, in mol/m3.  C0  dX  kC0 (1  X ) dt  (4.3)  where X is the fraction of conversion of the oxygen dependent reaction, C 0 is the initial concentration of O2 , in mol/m3. Heat release from oxygen dependent reactions can be expressed by  qox  C0VX r 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  dqox q  k  r H (C0V  ox )  k (C0V  r H  qox ) dt r H  (4.5)  dqox  A1e kt dt  (4.6)  Solving Eq. (4.5), we get  61  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.  Error  (  dqox dt   exp  dqox dt  )/( Eq.(06)  dqox dt  )  (4.7)  exp  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). k (T )  A0e E / RT  (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)  r H   A1M s kVC 0  (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  62  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.8x10 5 s), which corresponds to the heat release rate of the oxygen independent reactions. All tests were measured in two repeats.  0.4 50oC 40oC  0.35  Heat release rate (mW/g)  30oC 0.3 0.25  Oxygen dependent reactions dominate Oxygen independent reactions dominate  0.2 0.15 0.1 0.05 0  0  0.5  1  1.5  2  2.5  3  3.5  4  Time (s)  4.5 5  x 10  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  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 CO 2 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.  64  0.4 0.35 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%)  Heat release rate (mW/g)  0.3 0.25 0.2 0.15 0.1 0.05 0 -0.05 -0.1  0  0.5  1  1.5  Time (s)  2  2.5  3 5  x 10  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 65  (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.4 0.35 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%  Heat release rate (mW/g)  0.3 0.25 0.2 0.15 0.1 0.05 0 -0.05 -0.1  0  0.5  1  1.5  2  2.5  3  Time (s)  5  x 10  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  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.  67  Table 4.4 Reaction rate constant and activation energy of fresh samples with different moisture contents Temperature k (30 oC) *105, s-1  Sample 1  Sample 1  Sample 2  Sample 2  Sample 9  Sample 9  (Pellets B,  (Pellets B, M.C.  (Pellets B,  M.C. 4%)  4%)  M.C. 8%)  8%)  0%)  0%)  (Pellets B, M.C. (Pellets B, M.C. (Pellets B, M.C.  0.489±0.003  0.549±0.001  0.415±0.002  0.543±0.002  0.493±0.002  0.372±0.010  *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  k (40  oC)  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].  90  Activation energy (kJ/mol)  85 80 75 70 65 60 55 50 -1  0  1  2  3  4  5  6  7  8  x1 (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.  68  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 50 oC), 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. k (50o C )  2.063exp(0.014 x1 ) 105  (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.  69  2.2  -1  1.8  5  Reaction rate constant *10 (s )  2.0  1.6  o  k(30 C) o k(40 C) o k(50 C) Eq. (12)  1.4 1.2 1.0 0.8 0.6 0.4 0.2 -5  0  5  10  15  20  25  30  35  40  x1 (days) Figure 4.5 Comparison of experimental reaction rate constant data k(T) and Eq. (12) for group 2 pellets (stored at 23 oC). 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),. k (50o C )  1.894exp(0.255 x2 ) 105  (4.11)  where k(50oC) is the reaction rate constant at 50oC, and x2 is the days conditioned at 50oC before the test.  70  Group 2 Group 3 Eq. (10) Eq. (11)  2.0  5  -1  Reaction rate constant *10 (s )  2.2  1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.4  -5  0  5  10  15  20  25  30  35  40  x1 (days) Figure 4.6 Comparison of experimental reaction rate constant data and Eq. (4.10) and Eq.(4.11) for 50 oC 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 selfheating 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  k ( x1 , T )  2.883exp(0.014 x1   6801.78 ) 104 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.  71  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 50oC 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. Sample 1 Temperature  Sample 1  Sample 2  Sample 2  Sample 9  Sample 9  (Pellets B, M.C. (Pellets B, M.C. (Pellets B, M.C. (Pellets B, M.C. (Pellets B, M.C. (Pellets B, M.C. 4%)  4%)  8%)  8%)  0%)  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 (50oC), A1 value at 50oC, 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 A1 (50o C )  0.276exp(0.0256 x1 )  (R2=0.969) (4.13)  Group 3 A1 (50o C )  0.254exp(0.640 x2 )  (R2=0.982)(4.14)  where A1(50oC) is the heat release coefficient at 50oC, x1 is the days conditioned at 23oC before the test and x2 is the days conditioned at 50oC before the test.  72  0.30  Group 2 Group 3 Eq. (14) Eq. (15)  0.25  A1 (mW/g)  0.20  0.15  0.10  0.05  0.00 -5  0  5  10  15  20  25  30  35  40  x1 (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.  A1 ( x1 , T )  1.151exp(0.0256 x1   11637.80 ) 1015 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.  73  0.30  0.25  o  A1(30 C) o  A1(40 C)  A1 (mW/g)  0.20  o  A1(50 C)  0.15  0.10  0.05  0.00  -5  0  5  10  15  20  25  30  35  40  x1 (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 50oC, 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  74  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.  dq 11637.80 ( x1 , T )  1.151exp(0.0256 x1  ) 1015  dt T 6801.78 exp[2.883exp(0.014 x1  ) 104  t ]  0.0026 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 23oC 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.  75  0.4 Eq. (16) experimental data  0.35  Heat release rate (mW/g)  0.3  Fresh Pellets  0.25 10days 0.2 38 days 0.15 0.1 2 months 6 months  0.05  12 months  0  12 months -0.05 -0.1  0  0.5  1  1.5 2 Time (s)  2.5  3 5  x 10  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.15 Eq.(16) experimental data Fresh Pellets  0.1  Heat release rate (mW/g)  17 days 38 days 0.05  0 2 months one years old -0.05  -0.1  0  0.5  1  1.5 2 Time (s)  2.5  3 5  x 10  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.  76  0.15 Eq. (16) experimental data 0.1 Heat release rate (mW/g)  Fresh Pellets  0.05  0  2 months one years old  -0.05  -0.1  0  0.5  1  1.5  2 2.5 Time (s)  3  3.5  4  4.5 5  x 10  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.  77  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.  78  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.  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. 1  79  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 22 oC) Sample  C1  C2  Age  4oC, 3 months  4oC, 3 months  Head space  25%  0%  C3  C4  C5  C6  22oC, 3  22oC, 3  15oC, 2  15oC, 2  months  months  months  months  0%  25%  25%  0%  c. Sample group D for tests on age and head space effects (aged pellets prepared for tests at 40 oC)  Age  D1  D2  D5  D6  D7  D8  22oC , 3  22oC, 3  4oC, 3  4oC, 3  22oC, 2  15oC, 2  months  months  months  months  months  months  25%  0%  0%  25%  0%  0%  Head space  d. Sample group E for age and head space effects (aged pellets prepared for test carried out at 50 oC) 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  80  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:  fi   P(CiVg ) M wt CN 0 RTM PCNt  (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.  81  4.0 3.5  CO2 Concentration (%)  C13, 2L, 0% HS C14, 2L, 25% HS  A1, 45L, 0% HS B1, 45L, 25% HS  a  C1, 10L, 25% HS C2, 10L, 0% HS  3.0 2.5 2.0 1.5 1.0 0.5 0.0  CO2 Emission factors (g/kg)  b 0.03  0.02  0.01  0.00  0  10  20  30  40  50  60  70  Time (days)  Figure 5.1 Effects of reactor size and head space on (a) CO 2 concentration and (b) emission factors for fresh pellets at 22oC.  82  3.5  A1, 45L, 0% HS B1, 45L, 25% HS  CO Concentration (%)  3.0  C13, 2L, 0% HS C14, 2L, 25% HS  a  C1, 10L, 25% HS C2, 10L, 0% HS  2.5 2.0 1.5 1.0 0.5  CO Emission factors (g/kg)  0.0  b  0.020  0.015  0.010  0.005  0.000  0  10  20  30  40  50  60  70  Time (days)  Figure 5.2 The effects of reactor size and head space on CO concentration and emission factors for fresh pellets at 22 oC  83  18  Oxygen concentration (%)  C13,2L,0% HS C14,2L,25% HS  A1,45L,0% HS B1,45L,25% HS  16 14  C1,10L,25% HS C2,10L, 0% HS  12 10 8 6 4 2 0 -2  0  10  20  30  40  50  60  Time (days) Figure 5.3 The effects of reactor size and head space on O 2 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] Ci (t )  Ci , [1  exp(ki t )]  .  (5.2)  84  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 CO 2 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  fi (t )  fi , [1  exp(ki 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  85  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 CO 2 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 offgassing 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.  86  0.30  A1, 45L, 0% HS B1, 45L, 25% HS  CH4 Concentration (%)  0.25  C13, 2L, 0% HS C14, 2L, 25% HS  a  C1, 10L, 25% HS C2, 10L, 0% HS  0.20 0.15 0.10 0.05  CH4 Emission factors (g/kg)  0.00  b  0.0008 0.0006 0.0004 0.0002 0.0000  0  10  20  30  40  50  60  70  Time (days)  Figure 5.4 The effect of reactor size and head space on CH 4 concentration and emission factors for fresh pellets at 22oC  87  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 CO 2 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.  88  A1, 45L, 0% HS B1, 45L, 25% HS  1.8  H2 Concentration (%)  1.6  C13, 2L, 0% HS C14, 2L, 25% HS  C1, 10L, 25% HS C2, 10L, 0% HS  1.4 1.2 1.0 0.8 0.6 0.4 0.2  H2 Emission factors (g/kg)  0.0 0.0008 0.0007 0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0.0000  0  10  20  30  40  50  60  70  Time (Days)  Figure 5.5 The effect of reactor size and head space on H 2 concentration and emission factor for fresh pellets at 22oC  89  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 nonoxygen 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 curvefitting 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.  90  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 offgassing 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 4 oC 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 22oC. The emission kinetic parameters were fitted using the least square curve-fitting method following Eq. (5.3) and listed in Table 5.4.  91  CO2 Emission factor (g/kg(pellets))  0.030  C1 C2  0.025  C4 C3  C6 25% HS C5 0% HS  0.020  0.015  0.010  0.005  0.000  -0.005  0  10  20  30  40  50  60  70  Time (days) 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.  92  22  Oxygen concentration (%)  C4 C3  C1 C2  20  C6 25% HS C5 0% HS  18 16 14 12 10 8 6 4 2 0 -2  0  10  20  30  40  50  60  70  Time (days)  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 (22 oC) (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  93  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.  94  CH4 Emission factor (g/kg(pellets))  0.0007  C1 C2  0.0006  C4 C3  C6 25% HS C5 0% HS  0.0005 0.0004 0.0003 0.0002 0.0001 0.0000 0  10  20  30  40  50  60  70  Time (days) 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.  H2 Emission factor (g/kg(pellets) )  0.0008  C1 C2  0.0007  C4 C3  C6 25% HS C5 0% HS  0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0.0000 0  10  20  30  40  50  60  70  Time (days) 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.  95  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.  96  0.09  CO2 Emission factor (g/kg(pellets) )  o  0.08  o  o  22 C test  40 C test  50 C test  C2 C3  D5 D7  E2 E3  0.07  Fresh pellets Aged pellets  0.06 0.05 0.04 0.03 0.02 0.01 0.00 0  10  20  30  40  50  60  70  Time (days) Figure 5.11 Effect of temperature on CO2 emission factor for fresh and aged pellets  0.020 o  CO Emission factor (g/kg(pellets) )  22 C test  o  C2 C3  40 C test D5 D7  10  20  o  50 C test E2 E3  Fresh pellets Aged pellets  0.015  0.010  0.005  0.000 0  30  40  50  60  70  Time (days) Figure 5.12 Effect of temperature on CO emission factor for fresh and aged pellets  97  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 40oC 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 40 oC 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).  25 o  o  Oxygen concentration (%)  22 C test 40 C test 50oC test E2 Fresh pellets E3 Aged pellets  D5 D7  C2 C3  20  15  10  5  0 0  10  20  30  40  50  60  70  Time (days) Figure 5.13 Effect of temperature on O2 depletion for fresh and aged pellets  98  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.  99  CH4 Emission factor (g/kg(pellets) )  o  22 C test C2 C3  0.0016 0.0014  o  40 C test 50oC test E2 Fresh pellets D5 E3 Aged pellets D7  0.0012 0.0010 0.0008 0.0006 0.0004 0.0002 0.0000 0  10  20  30  40  50  60  70  Time (days) Figure 5.14 Effect of temperature on CH4 emission factor for fresh and aged pellets.  H2 Emission factor (g/kg(pellets) )  Temperature is added in the legend in the following figures according to dr. Bi's comment  o  o  o  0.0014  22 C test  40 C test  50 C test  0.0012  C2 C3  D5 D7  E2 E3  Fresh pellets Aged pellets  0.0010 0.0008 0.0006 0.0004 0.0002 0.0000 0  10  20  30  40  50  60  70  Time (days) Figure 5.15 Effect of temperature on H2 emission factor for fresh and aged pellets  100  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.  101  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.  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.. 1  102  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 84 oC 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 selfheating. 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 125 oC 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, selfheating 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  103  reactor is still lower than the oven temperature, which means that self-heating at this test segment is not sufficient to cause thermal runway.  SI11-19-2.xls  SI11-18-2.xls 110  100  90 105 Toven  o C)  o C)  80 T1 100  T2 T3  Temperature, (  Temperature, (  70  60 Toven T1 50  T2 T3  40  T4 95  90  T4 85  30  80  20 0  200  400  600  800  1000  1200  1400  0  1600  500  1000  1500  Time, (min)  Time, (min)  b  a  SI11-21-2.xls  SI11-20-2.xls 134  125  132  120  o C)  115  Temperature, (  Temperature, (  o C)  130  Toven 110  T1 T2 T3  105  T4  128  126  124  Toven T1  122  T2 T3  120  T4  100 118  116  95 0  500  1000  1500  2000  2500  0  3000  500  1000  1500  2000  2500  d  c SI11-24-2.xls  SI11-26-2.xls  240  250 Toven  Toven  T1  T1 220  T2  T2  200  T3  T3  T4  T4  o C)  o C)  3000  Time, (min)  Time, (min)  200  Temperature, (  Temperature, (  150  180  160  100  50 140  120  0 0  500  1000  1500  2000  Time, (min)  e  2500  3000  3500  0  1000  2000  3000  4000  5000  6000  Time, (min)  f  Figure 6.1 A set of measured temperature profiles from the self-heating reactor.  104  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. .  105  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.  106  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.  107  Table 6.1 The steady state center temperature for 3 sizes containers at different oven temperatures (Pellet E).  Oven  Large reactor  temperature  FS  100oC  SS  130oC  SS  140oC  SS  CP  Medium reactor  Max T  FS  100 oC yes yes  Max T  FS  SS  100 oC  SS  100 oC  132 oC  SS  130 oC  SS  130 oC  144 oC  SS  yes  141 oC  SS  140 oC  SS  yes  148 oC 150 oC  145 oC  CP  Small reactor CP  Max T  150oC  SS  162 oC  SS  153 oC  SS  160oC  TR  196 oC  SS  169 oC  SS  170oC  TR  -  SS  178 oC  SS  179 oC  180 oC  TR  -  TR  195 oC  SS  188 oC  191 oC  TR  -  TR  -  TR  200 oC  yes  162 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 130 oC; and after the crossing point, a stationary sate is reached with a maximum temperature of 132 oC 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 160 oC. 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.  108  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.  146  135  144  130  Temperature, (oC)  Temperature, (oC)  142  125  Toven L-T1 L-T2 M-T1 M-T2 S-T1 S-T2  120  115  110  138 Toven L-T1 L-T2 M-T1 M-T2 S-T1 S-T2  136 134 132  105  100  140  130 128  0  100  200  300  400  500  600  700  800  900  1000  0  200  400  165  200  160  190  155  150 Toven L-T1 L-T2 M-T1 M-T2 S-T1 S-T2  145  140  0  500  1000  1500  800  1000  1200  b. oven temperature 140oC  Temperature, (oC)  Temperature, (oC)  a. oven temperature 130oC  135  600  Time, (min)  Time, (min)  2000  Time, (min)  c. oven temperature 150 oC  180  170  Toven L-T1 L-T2 M-T1 M-T2 S-T1 S-T2  160  150  2500  140  0  200  400  600  800  1000  1200  1400  1600  Time, (min)  d. oven temperature 160 oC  Figure 6.3 Typical temperature profiles of three different sizes of reactors with increasing the oven temperature.  109  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  CP  T  2T   2    r h  A exp(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     2T    r h  A exp( E / RT ) x2  (6.2)  By introducing two dimensionless factors,   ( E / RT02 )(T  T0 ) ,   r / rc , Eq. (6.2) is converted to a dimensionless form,   2   e  2  (6.3)  in which the dimensionless parameter δ is defined as  110  E r hA 2  RT0  re RT02  E  (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  eff  (0.219  0.01M )  (1   )    0.027  (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  ln   2ln(T0 / r )  ln(  E r hA )  E / RT0 R .  (6.6)  The values of E and ∆rhA are obtained from a linear regression by plotting ln( C )  2ln(T0 / r ) as y coordinate and 1/ T0 as x coordinates. Figure 6.4 shows the linear regression of pellets A and Pellets E using FK method.  111  19.5  19.5  19.0  Pellet E y= -9849x+40.40 2 R =0.993  18.5  ln(C)+2ln(T/r)  ln(C)+2ln(T/r)  19.0  18.0  17.5  18.5  Pellet A 18.0  17.5  0.00217  0.00224  0.00231  0.00238  y= -9647x+40.43 2 R =0.995 0.00217  0.00224  1/T0  a. Pellet E Figure 6.4 Linear regression of  0.00231  0.00238  1/T0  b. Pellet A  ln( C )  2ln(T0 / r ) vs 1/ T0 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   CP  T t      r h  A exp( E / RTP )  (6.7)  P  where p indicates the crossing point; TP is the crossing point temperature, K; T / t P is the temperature increase rate at the crossing point, K/s; The value TP and T / t P are experimentally determined. By applying a natural logarithm, Eq. (6.7) becomes  112  ln(   hA T E )  ln( r )  t P CP RTP  (6.8)  Following equation (6.8), activation energy E is obtained from a linear regression by plotting ln(dT / dt ) P vs 1/ TP , 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  CP,1 = (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.  -5.0  Pellet E y= -8123x+11.77 2 R =0.537  ln(dT/dt)/ln(K/s)  -5.5 -6.0 -6.5 -7.0 -7.5 -8.0 -8.5 0.0021  0.0022  0.0023  0.0024  0.0025  1/Tp (1/K)  a. Pellets E Figure 6.5 Linear regression of ln(dT / dt ) P vs  b. Pellets A  1/ TP  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 (R 2), 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  113  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  * single step curve fitting, temperatures range from 100 to 180  oC  (oven temperature 170  0.95 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). 114  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 170 oC. 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×106, 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  115  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)*106 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.  116  300  300 280 260  Temperature, (oC)  Temperature, (oC)  250  200  150  Pellet A Pellet A Pellet C Pellet C Numerical Numerical  100  240 220 200 180 160  Pellet B Pellet B Numerical Numerical  140 120  50  100 0  0  100  200  300  400  500  600  700  800  900  80  1000  0  100  200  300  Time, (min)  400  500  600  700  Time, (min)  a. Pellets A and C in medium reactor at oven temperature  b. Pellet B in medium reactor at oven temperature 170 oC  170 oC  (initial temperature 100 oC)  300  300 280 Pellet E Pellet E Numerical Numerical  260  Temperature, (oC)  Temperature, (oC)  250  200  150 Pellet E Pellet E Numerical Numerical  100  240 220 200 180 160 140  50  120 0  0  100  200  300  400  500  600  700  800  900  1000  100  0  Time, (min)  1000  2000  3000  4000  5000  Time, (min)  c. Pellet E in large reactor at oven temperature 170 oC  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)  117  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  ln   2ln(T0 / r )  40.4  9748 / T0  .(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 twodimensional 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  118  z  L=2R  B.C Twall=Ta, at z=L; Twall=Ta, at z=0; Twall=Ta, at r=R; əT/ər=0, at r=0; R  r  Equation:  CP  T 1  T  2T E  [ (r )  2 ]     r h  A exp( ) t r r r z RT  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.  250  o  Critical Temperature ( C)  200  150  Wood fibre FK method prediction Experiment data 2D Numerical prediction  100  50  0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5  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;  119  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.  120  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 200 oC, 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 selfheating rate constant (∆rhA) of (4.22±2.5)*106 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.  121  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).  122  Figure 7.1 Illustration of the cylindrical silo and two dimensional axi-symmetric model  123  We chose one silo to simulate the self-heating process. Assuming the system is axisymmetric, a two dimensional axi-symmetric model is developed in this simulation. Axisymmetric 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 Qwall  hWall A(TW ,S  Tinf )  (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TP  W TW .  (7.2)  124  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. 4 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)  H  B. 1 B. 3  LW  R=D/2  B. 2 LW Case B : with airflow Boundary Condition for Air: B. 4 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  B. 1 Boundary Condition for Soild: B. 3 B.1 Axial symmetry at r=0. B.2 Thermal insulation, B.3 Newton cooling law Eq.(7.1) B.4 Thermal insulation,  B. 2  Figure 7.2 Diagram of computation domain and boundary condition for both cases with and without airflow  125  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)  1CP,1 where  1  T1 1  T1  2T1  1 ( (r )  2 )  QG (T1 )  QH t r r r z  (7.3)  is the bulk density of the pellets, kg/m3; C P,1 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/m3; r and z are the cylindrical coordinates, m; t is the time, s. For gas phase (Air)  2CP,2 (  Energy  Momentum  T2 T 1  T2  2T2  z 2 )  2 ( (r )  2 )  QH t z r r r z  2 z   1   z  2 z P (  z z )  2 ( (r )  2 )    2 g z  t z  r r r z z  2  (7.5)  2 ( 2z )  0 t z  Continuity  where  (7.4)  (7.6)  is the bulk density of the air, kg/m3; C P,2 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.  126  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)  dQ=hloc(aSdz)(T-T2) Convection and Conduction (chcc) (Gas phase)  Input {‘R’,T_init’,’T_i nf’’,’v0'}  Output {‘T’,’T2',’v’,’p’}  {‘T2’,’v’,’p’} Incompressible Navier-Stokes(chns) (Gas phase)  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 C P,1 ) have been  127  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.01M )  (1   )    0.027  (7.7)  CP,1 = (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.  11637.80 ) T exp[2.883exp(0.014 x1  6801.78 / T ) 1e4  t )]  0.0026   QG (T )   1.15e15exp(0.0256 x1   (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  QG (T )   1.15e15exp(0.0256 x1   11637.80 )  0.0026 T  (7.10)  where QG is the heat generated from wood pellets, W/m3. ρ 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. 128  In Chapter 5, the heat generation rate of wood pellets at high temperatures (100 oC 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. QG (T )    4.22 109  exp[78700 / (8.314T )] when 323K<T<473K  (7.11)  where, QG is the heat generated from wood pellets, W/m3. ρ 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.  Pellets  Twall  H  Tinf Pellets D Twall Tinf  Free Convection  Force convection Wind condition  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 129  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]  Nu free  C (Pr, shape)(Gr Pr)1/4  (7.12)  where, C is a function of Pr and the shape of heat transfer surface and can be expressed as  C  C1 (shape)C2 (Pr)  (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   g  H 3 (T0  T1 ) CP  (Gr Pr)  ( )( ) 2   (7.14)  where,  is the gas density at temperature T . β is the coefficient of volume expansion which is defined as    1 V 1  ( )P   ( )P . V T  T  Inserting the air properties, the heat transfer coefficient can be expressed as 2  hWall ,1  Nu free   gCP  3 1/4 k T 1/4  C1C2C3[( ) ]C4 [( ) ] H  TH  (7.15)  130  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.  Heat transfer coefficients (W/(m2 K))  4  3.5  3  2.5  2  1.5  H=1m H=5m H=10m H=20m  1  0.5 290  300  310  320  330  340  350  360  370  380  Wall surface temperature (K)  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 131  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  Nuconv  Nu free  0.6 Re1/2 Pr1/3  (7.16)  Heat transfer coefficients (W/(m2 K))  8  7  D=1m D=5m D=10m D=20m  6  5  4  3  2  1 0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2  Wind velocity outside of silo (m/s)  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 132  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.  dz  hloc is defined for a representative volume Sdz of the packed bed following [131]  Wood Pellets (T1)  dQ  Air flow (T2)  dQH  hloc ( AS Sdz)(T1  T2 )  (7.17)  where QH is the convective heat flux from wood pellets to airflow, W/m3. 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 jH  2.19 Re2/3  0.78Re0.381  (7.18)  Where jH, the Chilton-Colburn factor, and Reynolds number are defined as  jH   hloc CP  2/3 ( ) CPG0 k  (7.19)  133  Re   DS G0 6G0  (1   )  AS   (7.20)  Here, G0   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. A S 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  Convective heat transfer coefficients (W/(m2K))  hloc  0.663(G0CP )1/3 ( AS k )2/3  0.394G00.619 ( AS )0.381  0.286 (CP k 2 )1/3  (7.21)  140  120  100  80  T=273K T=373K T=473K  60  40  20  0  0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  Velocity of the forced ventilation inside silo (m/s)  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.  134  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,  p  a  b 2 L  (7.23)  (1   )2  where a  150 3 2 2  d  b  1.75  1   2 d 2  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.  135  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.22 20% 40% 60% 80%  Inside airflow velocity (m/s)  0.2 0.18  effciency effciency effciency effciency  0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02  0  2  4  6  8  10  12  14  16  18  20  Height of the bulk pellets (m)  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.  136  a.First stage, T init=300K  Difference  200 150 100 50 0 60s  360s  3600s  2h 4h Timestep b.Second stage, T init=300K  6h  12h  Difference  0.8 0.6 0.4 0.2 0 60s  360s  3600s  360s  3600s  2h 4h 6h Timestep c.Third+Fourth stage, T init=300K  12h  24h  12h  24h  Difference  30 20 10 0 60s  2h 4h Timestep  6h  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  Residue   t  end r  R z  H     (T  T ) t 1 r  0 z  0  1  0  2  (7.24)  137  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.  138  500 T init=310K  Center temperature (K)  480  T init=315K  460  T init=320K  440 420 400 380 360 340 320 300 0  10  20  30  40  50  60  70  80  Time (days)  Ae Time=3days  Be Time=10days 500  500  6  6  4  400  2 0  0  450  z,meter  z,meter  450 4  350  2  300  0  5 r,meter Ce Time=20days  400 350 300 0  5  r,meter De Time=25.5days 500  500  6  6  4  400  2 0  0  5 r,meter  450  z,meter  z,meter  450 4  350  2  300  0  400 350 300 0  5 r,meter  Figure 7.10 Four stages of self-heating.  139  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.  140  a. Outside wall temperature of different thickness and material  Temperature,K  305 304  303 Lc=0.025m, Concrete Lc=0.025m, Steel Lc=0.1m, Concrete Lc=0.1m, Steel Lc=0.15m, Concrete Lc=0.15m, Steel  302  301  300  0  10  20  30  40  50  60  70  Time, days b. Inside wall temperature of different thickness and material 303  Temperature,K  302.5 302 301.5  Lc=0.025m, Concrete Lc=0.025m, Steel Lc=0.1m, Concrete Lc=0.1m, Steel Lc=0.15m, Concrete Lc=0.15m, Steel  301 300.5 300  0  10  20  30  40  50  60  70  Time, days  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 141  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.  355  r=0m r=8m r=9m  350  Temperature (K)  345 340 335 330 325 320 315 310  0  10  20  30  40  50  60  70  Time (days)  Figure 7.12 The temperature profile of steel wall at different wind speeds v=(0,1,2,3 ) m/s  314 313.5  Temperature (K)  313  v=0 m/s  312.5 v= 1 m/s v= 2 m/s  312 311.5  v= 3 m/s  311 310.5  Inside wall Outside wall  310 309.5  0  10  20  30  40  50  60  70  Time (days)  Figure 7.13 The temperature development of inner and outer walls atdifferent wind speeds  142  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 selfheating 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.  Critical ambient temperature (K)  380  360 Spontaneous Combustion 340  320 Safe Region 300 0  2  4  6  8  10  Radius (m)  Figure 7.14 Relationship between critical ambient temperature and radius of the silo  143  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  z  z 2 1   z  2 z P (  z ) ( (r )  2 )    2 g z   2 g  (T  T )  t z  r r r z z  (7.25)  where 2 is the bulk density of the air, kg/m3; υ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 (   P  2 g z ) and buoyancy force (  2 g  (T  T ) ). The buoyancy force caused by z  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.  144  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.  145  results 314K  results 319K  400  400 v=0.0001 m/s v=0.001 m/s v=0.002 m/s  390  380  A2  370  Temperature (K)  Temperature (K)  380  360  350  340  A1  370  360  350  340  330  330  320  320  310  0  10  20  30  40  50  v=0.001 m/s v=0.003 m/s v=0.004 m/s  390  60  70  80  310  0  10  20  Time (days)  60  70  80  results 335K 430  A3  v=0.02 m/s v=0.04 m/s v=0.055 m/s  420 410  v=0.19 m/s v=0.24 m/s v=0.29 m/s  420 410  400  Temperature (K)  Temperature (K)  50  b. Ambient temperature 319 K  results 329K  390 380 370 360 350  400 390 380 370 360 350  340  340  A4  330 320  40  Time (days)  a. Ambient temperature 314 K  430  30  330 0  2  4  6  8  Time (days)  c. Ambient temperature 329 K  10  320  0  2  4  6  8  10  Time (days)  d. Ambient temperature 335 K  Figure 7.15 Simulation results of airflow effect at different ambient temperatures and airflow velocities  146  A1 T=314k, v=0.0001m/s Time=10days  A2 T=314k, v=0.0001m/s Time=62days  420 10  400 380  5  360  z (m)  z (m)  10  420 400 380 5  360  340 0  340  320 0  5 10 r (m) A3 T=329k, v=0.02m/s Time=2days  0  320 0  5 10 r (m) A4 T=329k, v=0.055m/s Time=3days 420 10  400 380  5  360  z (m)  z (m)  10  420 400 380 5  360  340 0  320 0  5 r (m)  10  340 0  320 0  5 r (m)  10  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.  147  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).  Critical ambient temperature (K)  350  Spontaneous Combustion  340  330  Safe Region  320  310 0.0  0.5  1.0  1.5  2.0  2.5  3.0  3.5  Airflow velocity (m/s)  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.  148  380  400 v=0.001 m/s v=0.002 m/s v=0.003 m/s  v=0.004 m/s v=0.005 m/s v=0.006 m/s  370  Center Temperature (K)  Center Temperature (K)  390 380 370 360 350 340 330  360  350  340  330  320  320 310  310  0  10  20  30  40  50  60  0  5  10  Time (days)  15  20  25  30  Time (days)  a. Ambient temperature 314K, airflow every other b. Ambient temperature 319K, airflow every other 24 hours 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 o  C), 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.  149  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.  150  60 55  L11 L10 L9 L8 L7 L6 L5 L4 L3 L2 L1  50  Temperature, oC  45 40 35 30 25 20 15 10 0  5  10  15  20  25  days from 10th June, 2009  Figure 7.19 Temperatures of cable C1 inside a industrial silo (80% load)  335  13meter 8.6meter 6.4meter 4.2 meter Ground  330  325  Center temperature(K)  Center temperature(K)  335  320  315  310  13meter 8.6meter 6.4meter 4.2 meter Ground  330  325  320  315  0  5  10  15  20  Time (days)  a.V=0.06 m/s  25  30  310  0  5  10  15  20  25  Time (days)  b.V=0.03 m/s  Figure 7.20 Simulation results for airflow velocities of 0.06 m/s and 0.03 m/s.  151  30  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 (53 oC). 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 (57 oC) 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 152  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.  level1  level3  30  30 Center T7 Half T6 Wall T5  Center T17 Half T16 Wall T15  28  Temperature (oC)  Temperature (oC)  28 26 24 22 20 18  26 24 22 20 18  0  100  200  300  400  Time (hour)  a. Temperature profile of level 1 (lower)  500  0  100  200  300  400  500  Time (hour)  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).  153  30  30 Center Half Wall  26  24  22  26  24  22  20  20  18  18 0  100  200  300  400  Center Half Wall  28  Temperature (oC)  Temperature (oC)  28  500  0  100  Time (hours)  a. Temperature profile of level 1 (lower) Tambient =  21oC;  Tinitial  200  300  400  500  Time (hours)  b. Temperature profiles of level 3 (upper)  =18.5oC  Tambient = 22oC; Tinitial =20.5oC,  Figure 7.22 Temperature profiles from pilot size reactor simulation results with unlimited oxygen supply  30  30 Center Half Wall  26  24  22  26  24  22  20  20  18  18 0  100  200  300  400  Time (hours)  Center Half Wall  28  Temperature (oC)  Temperature (oC)  28  500  0  100  200  300  400  500  Time (hours)  a. Temperature profile of level 1 (lower)  b. Temperature profiles of level 3 (upper)  Tambient = 21oC; Tinitial =18.5oC  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), 154  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 90 oC (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.  155  SI11-18-2.xls  SI11-19-2.xls  100  110  90  105  Temperature (oC)  Temperature (oC)  80 70 60 50  Toven T1 experimental results T3 (center) experimental results T1 simulation results T3 (center) simulation results  40 30 20  0  200  400  600  800  1000  1200  1400  100  95  90 Toven T1 experimental results T3 (center) experimental results T1 simulation results T3 (center) simulation results  85  80  75  1600  0  500  Time (min)  1000  1500  Time (min)  a. Oven temperature 90oC  b. Oven temperature 105oC SI11-21-2.xls  SI11-20-2.xls 138  125  136  120  Temperature (oC)  Temperature (oC)  134  115  110  Toven T1 experimental results T3 (center) experimental results T1 simulation results T3 (center) simulation results  105  132 130 Toven T1 experimental results T3 (center) experimental results T1 simulation results T3 (center) simulation results  128 126 124  100 122  95  120  0  500  1000  1500  2000  2500  3000  0  500  1000  1500  2000  2500  3000  Time (min)  Time (min)  d. Oven temperature 132oC  c. Oven temperature 125oC SI11-24-2.xls 260  Temperature (oC)  240  Toven T1 experimental results T3 (center) experimental results T1 simulation results T3 (center) simulation results  220  200  180  160  140  120  0  500  1000  1500  2000  2500  3000  3500  Time (min)  e. Oven temperature 143oC Figure 7.24 Comparison between experimental results and simulation results in lab scale reactor  156  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.  157  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 twodimensional 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.  158  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 selfheating 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)*106 kJ/(kg s) were obtained for four white pellets samples from different pellet producers in British Columbia.  159  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. 160  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.  161  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. 162  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.  163  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-businessin-norwich/.  4.  Firefighters battle huge biomass fire at Port of Tyne 2011. [Acessed 27th Apr. 2012]. http://www.forestbusinessnetwork.com/10442/firefighters-battle-huge-biomassfire-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-plantno-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/articleFirefighters-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.  164  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. 19601966. 165  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. 789796.  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 selfheating processes. Annals of Microbiology, 2003. 53(4): p. 349-410.  166  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.  167  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 threeparallel-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. 225239.  73.  Cashdollar, K.L., Coal dust explosibility. Journal of Loss Prevention in the Process Industries, 1996. 9(1): p. 65-76. 168  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. 169  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-termreaction 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.  170  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-cementbased 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 twophase porous geomaterials. International Journal of Heat and Mass Transfer, 2009. 52(34): 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. 171  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.  172  Appendix A Model parameters A.1. Material properties Wood pellets  Properties  Expressions / Values  Thermal conductivity,  1  (0.219  0.01M )  (1   )    0.027  λ, W/(mK) Specific heat capacity,  CP,1 = (1- M /100)(0.00546T - 0.524) + 0.042M  CP, J/(kgK) Density, kg/m3 Self-heating rate, Q=Q1+Q2, J/(kgK)  650-700  Q1(T )   1.15e15exp(0.0256 x1   11637.80 )  0.0026  T  With enough oxygen  Q2 (T )    4.22 109  exp[78700 / (8.314T )] Self-heating rate, Q=Q1+Q2, J/(kgK) With enough oxygen  11637.80 ) T exp[2.883exp(0.014 x1  6801.78 / T ) 1e4  t )]  0.0026   Q1 (T )   1.15e15exp(0.0256 x1   Q2 (T )    4.22 109  exp[78700 / (8.314T )] Average particle diameter, DS,  0.0061  m Average particle length,  0.010  LS,m ∅  0.42  a  477  173  Air  Properties  Expressions  Specific heat  CP  0.0769T  1076.9  1.1e3  Density, kg/m3    P  28.8e3 / (8.314T )  1.22  Thermal    100.8616log10( abs (T ))3.7142  0.0251    7.887 1012  T 2  4.427 108  T  5.204 106  1.7e-5  capacity, CP, J/(kgK)  conductivity, 𝛌, W/(mK) Dynamic Viscosity, 𝛍, Pa·s  Steel and Concrete  Properties  Steel  Concrete  Thermal conductivity,  44.5  1.8  475  750  7850  2300  𝛌, W/(mK) Specific heat capacity, CP, J/(kgK) Density, kg/m3  174  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)  175  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.  176  Figure C.1 Diagram of the pilot size reactor and the locations of the thermocouples[133].  177  level0left  level1left 30 Center T2 Half T3 Wall T4  28 26  Temperature (oC)  Temperature (oC)  30  24 22 20 18  Center T7 Half T8 Wall T9  28 26 24 22 20 18  0  100  200  300  400  500  0  100  Time (hour) level2left  400  500  30 Center T12 Half T13 Wall T14  28 26  Temperature (oC)  Temperature (oC)  300  level3left  30  24 22 20 18  Center T17 Half T18 Wall T19  28 26 24 22 20 18  0  100  200  300  400  500  0  100  Time (hour)  200  300  400  500  Time (hour)  level4left  level5left  30  30 Center T22 Half T23 Wall T24  28 26 24 22 20 18  Temperature (oC)  Temperature (oC)  200  Time (hour)  Center T27 Half T28 Wall T29  28 26 24 22 20 18  0  100  200  300  Time (hour)  400  500  0  100  200  300  400  Time (hour)  Figure C.2 Temperature profile from experimental results of lab scale experiments (Left side thermocouples) 178  500  level0right  level1right 30 Center T2 Half T1 Wall T0  28 26  Temperature (oC)  Temperature (oC)  30  24 22 20 18  Center T7 Half T6 Wall T5  28 26 24 22 20 18  0  100  200  300  400  500  0  100  Time (hour) level2right  400  500  30 Center T12 Half T11 Wall T10  28 26  Temperature (oC)  Temperature (oC)  300  level3right  30  24 22 20 18  Center T17 Half T16 Wall T15  28 26 24 22 20 18  0  100  200  300  400  500  0  100  Time (hour)  200  300  400  500  Time (hour)  level4right  level5right  30  30 Center T22 Half T21 Wall T20  28 26  Temperature (oC)  Temperature (oC)  200  Time (hour)  24 22 20 18  Center T27 Half T26 Wall T25  28 26 24 22 20 18  0  100  200  300  Time (hour)  400  500  0  100  200  300  400  Time (hour)  Figure C.3 Temperature profile from experimental results of lab scale experiments (Right side thermocouples) 179  500  Figure C.4 O2 concentration inside of the silo during the experiments. (Reprinted from Yazdanpanah et. al. [132] with permission from Yazdanpanah F.)  180  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')];  181  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;  182  %------------------------ 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'};  183  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};  184  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};  185  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  186  % 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)]';  187  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.8e3/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';  188  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  189  Appendix E. Self-heating with free convection  During the self-heating simulation, it is assumed that the free convection effect on selfheating 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. 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)  H  B. 1 B. 3  LW  R=D/2  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.  B. 2 LW  Figure E. 1 Diagram of computational domain and boundary condition E.2 Governing equation For the solid phase (Wood Pellets) 190  1CP,1 where  1  T1 1  T1  2T1  1 ( (r )  2 )  QG (T1 )  QH t r r r z  (7.3)  is the bulk density of the pellets, kg/m3; C P,1 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/m3; r and z are the cylindrical coordinates, m; t is the time, s. For gas phase (Air)  2CP,2 (  Energy  Momentum  T2 T 1  T2  2T2  z 2 )  2 ( (r )  2 )  QH t z r r r z  2 z   1   z  2 z (  z z )  2 ( (r )  2 )   2 g  (T  T )  t z  r r r z  2  (E.1)  2 ( 2z )  0 t z  Continuity  where  (7.4)  (7.6)  is the bulk density of the air, kg/m3; C P,2 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  191  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.  Ae Time=3days  Be Time=10days 500  4  400  2 0  6  450  0  z,meter  z,meter  6  500 450  4  400  350  2  350  300  0  5 r,meter Ce Time=20days  300 0  5  r,meter De Time=23.2days 500 450  4  400  2 0  500 6  0  5 r,meter  z,meter  z,meter  6  450  4  400  350  2  350  300  0  300 0  5 r,meter  Figure E.2 Solid phase temperature profiles at different stages. 192  Figure E.3 Air velocity caused by the free convection at the 11th day of self-heating.  193  

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}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            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