MEASUREMENT OF THE THERMOCHEMICAL ENERGY DISCHARGE PROCESS IN REACTIVE PACKED BEDS OF MAGNESIUM-MANGANESE OXIDE By Keith King A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering - Doctor of Philosophy 2021 ABSTRACT MEASUREMENT OF THE THERMOCHEMICAL ENERGY DISCHARGE PROCESS IN REACTIVE PACKED BEDS OF MAGNESIUM-MANGANESE OXIDE By Keith King Decarbonization of the world’s grid-scale electricity infrastructure will require the implementation of a variety of energy storage technologies. One such technology is thermochemical energy storage with reactive packed beds of metal oxide materials that undergo cyclically stable reduction/oxidation (redox) reactions at high temperatures. These materials store energy via an endothermic reduction reaction that releases oxygen gas, then discharge energy by a coupled exothermic oxidation reaction and interstitial convective heat transfer process, with atmospheric air as both the reactant and heat transfer fluid. Magnesium-manganese oxide is a candidate redox metal oxide material that has demonstrated cyclical reactive stability between 1000-1500 °C for several compositional variants. In order to build an energy storage reactor based on a packed bed of these materials, it’s critical that the overall energy storage capacity, the optimal operating conditions, and the rate of energy release during the discharge step are quantified. The first half of this work develops a calorimetric method for determining the energy storage densities of magnesium-manganese oxides of five Mn/Mg compositional variants between 1000-1500 °C. Drop calorimetry measures the total enthalpy of the phases that exist at the drop temperature, from which the energy stored as sensible heat is computed. Acid solution calorimetry of samples rapidly-quenched in inert gas measures the standard enthalpy of formation at room temperature of the compounds that exist at high temperatures, from which the energy of the redox reaction is computed. The total energy density is then summed from these measurements. The further gains in energy density that can be achieved by lowering the partial pressure of oxygen during reduction are also examined. The Mn/Mg variant of 1/1 is found to have the maximum energy storage capacity when reducing under a low oxygen partial pressure and cycling between 1000-1500 °C. A bulk oxidation kinetic model is also presented as part of the energy discharge step. The second half of this work quantifies the convective heat transfer within a packed bed of magnesium-manganese oxide particles. Porous medium heat and mass transport theory for geometrically homogenous and isotropic packed beds is summarized. A collection of similar studies from the literature are selected for comparison with this work. Both a steady-state and a transient method are considered for measuring the interstitial convective heat transfer coefficients. The transient method is found to be more suitable, and a general particle-diameter based convection correlation is proposed for packed beds of spheres. Unique correlations are developed with the transient method for two beds of 1/1 Mn/Mg particles, one granular and one cylindrical. A parametric analysis of the effect of transient wall-to-fluid convective heat transfer on the accuracy of the measured interstitial heat transfer coefficients is performed, finding that a minimum bed radius-to-bed length (R/L) ratio between 1 and 2 should be chosen for any similar experiments and that a lower minimum Reynolds number requires a larger minimum R/L ratio in order to minimize the measurement error. Finally, a 1D coupled oxidation-convection model combines these measurements in a parametric study of the outlet behavior of a reactive 1/1 Mn/Mg packed bed for different inlet Reynolds number and length-to-radius ratios and determines the scenarios in which accurate interstitial convective heat transfer coefficient measurements are important for modeling and controlling the discharge step of a magnesium- manganese oxide energy storage reactor. ACKNOWLEDGEMENTS I’d like to thank all of these people for their roles in helping me reach this point. To Dr. Neil Wright and Dr. Wei Lai, for serving on my dissertation committee and for their valuable feedback on this work. To my advisors, Dr. James Klausner and Dr. Jörg Petrasch: for giving me this opportunity; for their wisdom, expertise, mentorship, and continual guidance over the course of my PhD; and most importantly, for the kindness and patience they both have shown and continue to show me as I have learned and grown as an engineer. To all of my lab friends and colleagues over the last four years: Dr. Kelvin Randhir, Dr. Amey Barde, Dr. Clement Roy, Dr. Nima Rahmatian, Joe Kerwin, Sarah Plant, Alessandro Bo, Phillip Schimmels, and Michael Hayes, for their continual advice and friendship, and especially to Kelvin, who has taught me so much and been a constant source of support since my very first day at Michigan State. And finally, to my friends and family, especially my parents, for their constant love, support, and for keeping me sane throughout this whole experience. I would not have made it through without all of you. iv TABLE OF CONTENTS LIST OF TABLES viii LIST OF FIGURES ix KEY TO SYMBOLS AND ABBREVIATIONS xii CHAPTER 1: INTRODUCTION TO THERMOCHEMICAL ENERGY STORAGE AND MAGNESIUM-MANGANESE OXIDE 1 1.1. Context and Motivation 1 1.2 Thermochemical Energy Storage 3 1.2.1. Fundamentals 3 1.2.2. Summary of Approaches 4 1.3. Redox Metal Oxides for High-Temperature TCES 6 1.3.1. Advantages and Review of Materials 6 1.3.2. The Manganese Redox System and Magnesium-Manganese Oxide 8 1.3.3. Thermophysical Properties of Materials for Reactor Design 13 1.4. Objectives and Organization 14 CHAPTER 2: THERMOCHEMICAL PROPERTIES OF MAGNESIUM-MANGANESE OXIDE 17 2.1. Introduction 17 2.2. Predictions from CALPHAD Models 18 2.3. Calorimetric Method for Measuring Energy Density 23 2.3.1. Introduction to Calorimetry with High-Temperature Metal Oxides 23 2.3.2. Calorimetric Approach 24 2.4. Experimental Setup and Procedure 28 2.4.1. The Drop Calorimeter 28 2.4.2 The Acid Solution Calorimeter 29 2.4.3. The Calorimetric Procedure 30 2.5. Calorimeter Calibration and Validation 33 2.5.1. The Drop Calorimeter 33 2.5.2. The Acid Solution Calorimeter 33 2.5.2.1. Addition of SnCl2 to the Calorimetry Process 33 2.5.2.2. Calibration & Validation 38 2.6. Measurements for Magnesium-Manganese Oxide 44 2.7. Optimizing for Maximum Volumetric Energy Density 50 2.7.1. Experimental Methods 52 2.7.1.1. Variation of Molar Composition 52 2.7.1.2. Lowering Oxygen Partial Pressure During Reduction 53 2.7.2. Results 57 2.8. Comparison of Energy Densities with CALPHAD Predictions 63 2.9. Bulk Oxidation Reaction Kinetics 66 v 2.9.1. Experimental Setup and Procedure 67 2.9.2. Kinetic Modeling and Data Extraction 68 2.9.3. Results 70 2.10. Summary 75 CHAPTER 3: OVERVIEW OF POROUS MEDIUM HEAT AND MASS TRANSFER THEORY 77 3.1. Volume Averaging Theory 78 3.2. Fluid Mechanics in Homogenous Porous Media with Single-Phase Flow 83 3.3. Conduction Heat Transfer in Homogeneous Porous Media 85 3.4. Convective Heat Transfer in Homogeneous Porous Media 86 3.5. Summary 89 CHAPTER 4: INTERSTITIAL CONVECTIVE HEAT TRANSFER IN PACKED BEDS OF MAGNESIUM-MANGANESE OXIDE 90 4.1. Review of the Experimental Literature 90 4.2. Steady State Measurement 91 4.3. Transient Measurement 93 4.3.1. Experimental Setup and Method 95 4.3.2. Thermal Effect of the Tube Wall and Volume-Averaged Model 102 4.3.3. Data Analysis: Numerical Method, Data Extraction, and Uncertainty 105 4.3.4. Results 111 4.3.5. Lumped Thermal Analysis of Particles 123 4.3.6. Scale Analysis of the Transient Wall Heat Transfer Effect 125 4.4. Summary 131 CHAPTER 5: IMPORTANCE OF INTERSTITIAL CONVECTION IN A REACTING PACKED BED OF MAGNESIUM-MANGANESE OXIDE 134 5.1. Coupled Model 134 5.2. Importance of Interstitial Convection 137 5.2.1. Inlet Reynolds Number 137 5.2.2. Aspect Ratio of Packed Bed 138 5.2.3. Discussion 139 CHAPTER 6: FUTURE WORK AND CONCLUSION 141 6.1.Future Work 141 6.2. Concluding Summary 142 APPENDICES 144 APPENDIX A: SUMMARY OF MATERIALS FOR THERMOCHEMICAL ENERGY STORAGE 145 APPENDIX B: SUPPLEMENTARY WORK FOR THE THERMODYNAMIC MEASUREMENTS OF MAGNESIUM-MANGANESE OXIDE 151 APPENDIX C: ATTEMPTS AT STEADY STATE INTERSTITIAL HEAT TRANSFER COEFFICIENT MEASUREMENTS 162 vi APPENDIX D: DERIVATION AND SCALING OF THE MODIFIED C-S MODEL 175 APPENDIX E: LOW-TEMPERATURE SPECIFIC HEAT MEASUREMENTS OF MAGNESIUM-MANGANESE OXIDE 186 APPENDIX F: INTERSTITIAL HEAT TRANSFER COEFFICIENT TRANSIENT MEASUREMENT DATASETS 187 APPENDIX G: STATISTICAL CHARACTERIZATION OF THE GEOMETRIC DISTRIBUTION OF THE GRANULAR PARTICLE BED 194 REFERENCES 200 vii LIST OF TABLES Table 1.1. TCES Chemical Systems 5 Table 2.1. Standard Reduction Potentials 34 Table 2.2. Free Energy of Dissolution of Me3+ Compounds 35 Table 2.3. Acid Solution Calorimetry Chemical Reactions 39 Table 2.4. Redox Condition Reference Numbers 52 Table 2.5. Oxidation Kinetic Parameters 70 Table 4.1. Instrumentation and Uncertainties 110 Table 4.2. Tube-to-Particle Diameter Ratios 119 Table B.1. Reactivity 159 Table B.2. Energy Density 160 Table B.3. Acid Solution Calorimeter Validation 160 Table B.4. Drop Calorimetry Results 160 Table B.5. Acid Solution Calorimetry Results 161 Table B.6. Mass Change Results 161 Table B.7. Energy Density Results 161 Table C.1. Instrumentation Error 165 viii LIST OF FIGURES Figure 1.1. EIA Grid-Scale Growth Projections 2 Figure 1.2. TCES System Schematic 4 Figure 1.3. Reactive Cycling Stability 12 Figure 2.1. Reactivity Predictions 21 Figure 2.2. Energy Density Predictions 22 Figure 2.3. Calorimetric Method Diagram 27 Figure 2.4. Drop Calorimeter 29 Figure 2.5. Acid Solution Calorimeter 30 Figure 2.6. SnCl2 Dissolution Enhancement, Predictions 35 Figure 2.7. SnCl2 Dissolution Enhancement, Test 37 Figure 2.8. Acid Solution Calorimeter Validation 43 Figure 2.9. XRD Data of Quenched Samples 45 Figure 2.10. Drop Calorimetry Results 47 Figure 2.11. Acid Solution Calorimetry Results 48 Figure 2.12. Excess Oxygen Content 49 Figure 2.13. Energy Density Results 50 Figure 2.14. Illustration of Reactivity for Varied PO2 54 Figure 2.15. Bulk Densities 56 Figure 2.16. Reactive Stability for New Mn/Mg Samples 57 Figure 2.17. Measured Reactivity for Varied PO2 59 Figure 2.18. Volumetric Energy Densities at Varied PO2 60 ix Figure 2.19. Comparison with CALPHAD, Reactivities 64 Figure 2.20. Comparison with CALPHAD, Energy Densities 65 Figure 2.21. Oxidation Kinetic Experiments, Setup 68 Figure 2.22. Oxidation Kinetic Experiments, Total Absorbed 71 Figure 2.23. Oxidation Kinetic Experiments, Data and Model Fits 72 Figure 3.1. Representative Elementary Volume 79 Figure 3.2. Intrinsic and Superficial Averages 81 Figure 4.1. Transient Experimental Setup, Test Section 98 Figure 4.2. Transient Experimental Setup 100 Figure 4.3. Transient Wall Heat Transfer Effect 103 Figure 4.4. Data and Simulated Best-Fit Curves 112 Figure 4.5. Increased Error at Lower Flow Rates 114 Figure 4.6. Results, Nusselt 116 Figure 4.7. Results, Chilton-Colburn 120 Figure 4.8 Results, Chilton-Colburn, Spheres 121 Figure 4.9. Results, Nusselt, 1/1 Mn/Mg Particles 122 Figure 4.10. Lumped Thermal Analysis 124 Figure 4.11. Dimensionless Error Analysis 129 Figure 5.1. Effect of Error, Reynolds 138 Figure 5.2. Effect of Error, L/R 139 Figure B.1. Dissolution Experiments, Known Materials 153 Figure B.2. Dissolution Experiments, Magnesium-Manganese Oxide 155 Figure C.1. Experimental Setup 168 x Figure C.2. Datasets 169 Figure C.3. Results, Nusselt 171 Figure C.4. Results, Chilton-Colburn 172 Figure E.1. DSC Measurements 186 Figure F.1. Transient Interstitial HTC Outlet Data and Best-Fit Simulations 188 Figure G.1. Image, with Scale 195 Figure G.2. Processed Image 196 Figure G.3. Histogram of Effective Spherical Diameter Distribution 197 Figure G.4. Box Plot of the Effective Spherical Diameter Distribution 198 xi KEY TO SYMBOLS AND ABBREVIATIONS Variables a Specific surface area (m-1) A,B,C Generic chemical reactant/product chemical formulae b Forcheimer friction factor (-) Bi Biot number (-) cp Specific heat (J kg-1 K-1) d Diameter (m) D Tube diameter (m) Ddzz Axial dispersion coefficient (m2 s-1) E Energy (kJ) Ė Time rate of energy change (W) Ea Activation energy (kJ mol-1) Eo Standard reduction potential (V) G Body forces (N) G Gibbs free energy (kJ mol-1) H Enthalpy (kJ mol-1, J kg-1 , kJ kg-1) h Heat transfer coefficient (W m-2 K-1) I Current (A) jH Chilton-Colburn factor (-) k Thermal conductivity (W m-1 K-1) K Permeability (m2) xii Ko Kinetic rate constant (s-1) L Length (m) m Mass (g, kg) m Slope ṁ Mass flow (kg s-1) m" Mass flux (kg m-2 s-1) M Molecular weight (g/mol, kg/mol) Me Generic metal elemental symbol n Normal vector to a surface (-) n Number of moles (mol) N Number of points (-) NRMSE Normalized root mean squared error (-) Nu Nusselt number (-) P Pressure (Pa, atm, bar) Ṗ Power (W) Pr Prandtl number (-) q” Heat flux (W m-2 K-1) Q̇ Heat flow (W) r Radial distance (m) r Ratio (-) ṙ Rate of species production (kg m-3 s-1) R Universal gas constant (J mol-1 K-1) xiii R Outer radius (m) Re Reynolds number (-) S Kinetic parameter (-) St Stanton number (-) T Temperature (°C, K) t Thickness (m) t Time (s) u Axial velocity (m s-1) u" Velocity vector (m s-1) v̇ Flow rate (SCCM, SLPM) V Voltage (V) V Volume (m3) V Volume of gas released per unit mass of solid (cm3 g-1) V̇ Rate of volume of gas released per unit mass of solid (cm3 g-1 min-1) x Mass fraction (-) x Ratio of magnesium/manganese (molMg molMn-1) y Excess oxygen atoms (molO mol-1MgxMnO1+x+y) z Axial distance (m) Greek Variables/Symbols 𝛼 Extent of reaction (-) 𝛼 Thermal diffusivity (m2 s-1) 𝛾 Mole fraction (-) xiv 𝛿 Extent of non-stoichiometry (-) Δ Change/difference in a value between two points 𝜖 Bulk porosity (-) 𝜖 Phase fraction (-) θ Dimensionless temperature (-) 𝜇 Uncertainty in a variable μ Dynamic viscosity (Pa-s) 𝜌 Density (kg m-3) 𝜏 Time constant (s) 𝜙 Generic property Superscripts and Subscripts * Signifying a dimensionless quantity ∞ Value at end time/equilibrium state 0 Initial value air Value specifically for air aspect Related to an aspect ratio Avg Average value bed Related to the packed bed scale CV Control volume D Darcy DC Direct current Diss Dissolution xv Drop Related to drop calorimetry eff Effective end Value at the end point in time of a dataset exp Experimental value f Formation f Fluid fit Curve fit value HTC Heat transfer coefficient i Index j Index in Inlet Int Interstitial Loss Related to heat loss mod Modeled/simulated value O2 Value related to oxygen gas out Outlet OX Oxidized material state p Particle RED Reduced material state Ref Reference value Rxn Reaction S Stored s Solid xvi sf Related to surface area interface between solid and fluid phases Sys System total Summed value w Wall Abbreviations DC Direct current DSC Differential scanning calorimetry HTC Heat transfer coefficient SCCM Standard cubic centimeters per minute SLPM Standard liters per minute TCES Thermochemical energy storage xvii CHAPTER 1: INTRODUCTION TO THERMOCHEMICAL ENERGY STORAGE AND MAGNESIUM-MANGANESE OXIDE 1.1. Context and Motivation The environmental and economic challenges presented by anthropogenic climate change and the accompanying need for decarbonization of the world’s energy economy are firmly established. Current global emissions targets under the 2015 Paris Climate Agreement have set many nations on the path to decarbonized or carbon-neutral economies by 2050, with steep goals as early as 2025-2030 in some developed nations [1]. As part of this agreement, the United States is obliged to reduce its green-house gas emissions by 26-28% of 2005 levels by 2025 [2]. Further, at the time of this writing, the current Biden administration has set targets for national net-zero emissions no later than 2050 [3]. As part of larger efforts towards meeting these goals, American grid-scale power generation infrastructure requires a rapid pivot from fossil fuels to carbon-free technologies, such as renewable-plus-storage systems. This is certainly a difficult task: in 2019, fossil fuels accounted for 62.6% of the roughly 4.13 trillion kWh of grid-scale electricity generated in the United States, whereas only 8.8% came from wind and solar photovoltaics (PV) [4] . However, wind and solar PV are currently experiencing swift growth in their share of total grid-scale generation capacity. Figure 1.1 from the United States Energy Information [5],[6] Administration (EIA) shows the projected growth in total installed generation capacity in the United States over the next three decades until 2050. 1 Figure 1.1. EIA Grid-Scale Growth Projections. EIA projections for growth of grid- scale electricity generation capacity in the United States until 2050, broken down by energy source [6]. Total installed renewable capacity is project to grow faster than all other sources, with new wind and solar PV installations making up almost all of the projected growth. If the EIA’s projections are correct, wind and solar PV will meet an increasingly larger portion of total grid-scale electricity demand over the next three decades. However, it is well- established that due to the intermittent nature of wind and solar energy, these technologies will require accompanying energy storage systems if they are to continue growing and become a foundational part (among other carbon-free generation technologies) of a decarbonized electric grid. Thus, opportunities exist for incorporating energy storage technologies into these new installations. The University of Michigan’s Center for Sustainable Systems website (consisting of information taken from the U.S. Dept. of Energy’s Global Energy Storage Database Projects 2020) details the total storage capacity and the percentage of new storage projects in 2020 across different technologies [7]. On evidence, the storage market is currently dominated by pumped hydropower at 94% of total installed storage capacity. However, because pumped hydro is limited in application by geography, there are opportunities for alternative storage technologies that are more suitable for regions that are more geographically diverse. 2 One such option is the implementation of thermal storage systems. Most established thermal storage systems are based on sensible heat storage in which some material such as water, an oil, or a molten salt (generally a mixture of nitrates and nitrides) is directly heated, taking advantage of the material’s thermal mass to store energy as heat [8]. Often coupled with concentrated solar power plants, these systems store relatively small amounts of energy per unit mass or per unit volume (i.e., energy density) of material, requiring both a large material mass and significant insulation to be economically practical [8]. Other existing thermal storage systems use the latent heat associated with a melting/solidifying phase change cycle as a storage mechanism to improve energy density [8]. Additional gains in energy density are possible by using reactive materials in a thermal storage system. This approach is defined as thermochemical energy storage. 1.2 Thermochemical Energy Storage 1.2.1. Fundamentals Thermochemical energy storage (TCES) uses a reversible chemical reaction to store energy in a two-step process. Energy input as heat drives an endothermic reaction of some reactive material in a “charging” step, creating two products that are isolated and stored separately until the energy is needed, at which point the products are re-combined in an exothermic reaction, releasing the stored energy in a “discharging” step [8]-[11]. The simplest form of this reaction is A↔B+C ΔHRxn > 0 (1.1) where A, B, and C are generic chemical formulas for the storage material reactant and products, respectively, and ΔHRxn is the enthalpy of the reaction. Adding this chemical reaction energy to the sensible & latent energy storage increases the overall energy density of the storage material. 3 Further discussion of the advantages of TCES compared with thermal and latent-heat storage systems can be found in [8]-[11]. An extremely simplified schematic of a TCES system can be seen in Figure 1.2. Figure 1.2. TCES System Schematic. Schematic of a simple theoretical thermochemical energy storage system [11]. The reactant A is decomposed to products B and C in an endothermic reaction driven by an energy source (solar energy in this schematic). Products B and C are re-combined in an exothermic reaction to drive electricity production when demand is present (symbolized by the power plant illustration). 1.2.2. Summary of Approaches Many chemical systems have been studied for TCES purposes, all of which consist of gas-gas, gas-liquid, or gas-solid reaction couples. These include ammonia, metal carbonates, metal hydrides, hydrocarbons (particularly methane), metal hydroxides, and sulfur oxides. Table 1.1 gives the fundamental storage reaction for each system and the accompanying enthalpy of reaction. A summary of the benefits and drawbacks of each system can be found in Appendix A. 4 Table 1.1. TCES Chemical Systems. Summary of the different chemical reactions under study for use in thermochemical energy storage systems, excluding reduction/oxidation metal oxides [8]-[11]. The recant name, basic stoichiometric reaction used to store energy, and the reaction enthalpies (per unit mass of the reactant) are shown. Chemical Basis Fundamental Reaction ΔHRxn (kJ kg-1) Ammonia 2NH3 (g) ↔ N2 (g)+3H2 (g) 3922 1 Hydrides MeH(s) ↔ Me (s)+ H2 (g) 1158-8397 2 Hydroxides Me(OH)2 (s) ↔ MeO (s)+H2 O (g) 1340-1790 CH4 (g)+H2 O (g) ↔ 3H2 (g)+CO (g) 1.54E4 Hydrocarbons CH4 (g)+CO2 (g) ↔ 2H2 (g)+2CO (g) 1.56E4 Carbonates MeCO3 (s) ↔ MeO(s)+CO2 (g) 500-1790 Sulfur 2SO3 (g) ↔ 2SO2 (g)+O2 (g) 1224 The chemical systems in Table 1.1 all require one or both of the following: (1) the need for separate reactors for the endothermic and exothermic steps, and (2) a gas storage and handling system to contain and transport the reactive gases. These requirements will contribute to the total cost of a TCES system based on one of these compounds. When the environmental health and safety concerns of each compound are also considered, the cost of building a TCES system around one of these storage materials could become prohibitively high [8]-[11]. A class of materials offering a solution to this challenge are metal oxides that undergo reversible reduction-oxidation (redox) reactions at high temperatures. The generalized reaction is given by Eq. 1.2 as y MeOx (s) ↔ MeOx-y (s) + 2 O2 (g) ΔHRxn > 0 (1.2) [8]-[11] 5 where Me is a generic metal atom, MeOx is the generic metal oxide in the oxidized phase, MeOx-y is the generic metal oxide in the reduced phase, x is the number of moles of oxygen atoms present in the oxidized phase, and y is the number of moles of oxygen atoms released during reduction. The heat applied to the metal oxide reduces the metal cations, allowing some oxygen to escape from the crystal lattice and recombine in a gaseous state. The reduced state then remains until the material is exposed to gaseous oxygen at lower temperatures and/or in higher concentrations, at which point oxygen atoms are absorbed into the lattice and heat is released in an oxidation reaction, returning the solid to its original crystal structure [8]-[11]. Materials that exhibit this reversible redox reaction are promising candidates for TCES applications and are briefly reviewed in the next section. 1.3. Redox Metal Oxides for High-Temperature TCES 1.3.1. Advantages and Review of Materials The biggest advantage of redox metal oxides is the ability to utilize ambient air as both heat transfer fluid and reactant. Consequently, these compounds operate as open systems, drawing and discharging gas to and from the atmosphere, thereby avoiding the need for storage tanks for potentially toxic reaction products [8]-[11]. The ideal material would exhibit all of the following characteristics [8]-[11]: 1. Inexpensive. 2. Reasonably abundant on Earth. 3. Minimal environmental concerns (non-toxic, non-flammable, non-explosive, non- corrosive, etc.). 4. Moderate-to-high reaction enthalpies and energy densities. 5. Suitably high (T > 700 °C) reaction temperatures. 6 6. Full reversibility and long-term cyclical stability with low sintering. 7. Stable mechanical and thermo-physical properties. 8. Minimal thermal hysteresis between reduction and oxidation temperatures. 9. No by-products or negative side-reactions. Compounds that have been studied include oxides of a range of elements from several sections of the periodic table, including lithium (alkali metal), barium, magnesium (alkaline-earth metals) chromium, copper, cobalt, iron, manganese, platinum, rhodium, vanadium (transition metals), lead (metal), antimony (semi-metal), and uranium (actinide metal). Of these, the barium, cobalt, copper, iron, and manganese oxides possess enough of the desired characteristics to warrant further investigation [8]-[11]. It should be noted that while barium is an alkaline-earth metal and the remainder (cobalt, iron, copper, manganese) are transition metals, all of their oxide compounds will be referred to as “redox metal oxides” for simplicity in the remainder of this work. The fundamental reaction couple, theoretical energy density, and some benefits and drawbacks of each can be found in Appendix A. A variety of mixed metal oxide compounds have also been studied as potential TCES materials. A mixed metal oxide compound combines another metal oxide with one of the pure redox metal oxides in an attempt to improve the TCES properties of the material. This has achieved varying degrees of success in enhancing energy density, increasing reaction kinetics, improving reversibility, and increasing cyclic stability by reducing sintering. In many cases, the mixed metal oxide systems that succeed in improving TCES properties often correspond to a loss in energy density relative to the pure compound [8]-[11]. Thus, these mixed systems come with trade-offs, and an optimal balance of energy density, reactive stability, reversibility, and sufficient reaction kinetics must be found. An extensive review of the work done studying both 7 pure and mixed redox materials can be found in [8]-[11] and the subsequent literature. In regards to its relevance to this work, further discussion of the manganese redox system will be given in the next section. Of the mixed metal oxides, it is worth mentioning the perovskite class of materials due to the fundamental difference in their reactive behavior. Perovskites are unique from the others because they do not change phase between two well-defined compositions during reduction /oxidation, but rather reduce gradually with increasing temperature once a threshold temperature has been crossed at a given oxygen partial pressure. Vacancies steadily form as the temperature is increased and oxygen is released non-stoichiometrically. Comparably, the oxidation reaction occurs steadily as the temperature is reduced and heat is gradually released, with the overall energy density proportional to the extent of oxidation [8],[9]. These materials have the generic chemical formula of ABO3 and exhibit the generalized reaction δ ABO3 (s) ↔ ABO3-δ (s)+ 2 O2 (g) ΔHRxn > 0 (1.3) [8] where A and B are different metal cations and 𝛿 is the non-stoichiometric extent of reduction. A range of perovskites have been examined with different combinations of lanthanum, strontium, cobalt, manganese, iron, barium, calcium, chromium, yttrium, aluminum, and titanium on the A and B sites in the lattice [8]-[9]. These materials are promising candidates for further exploration and tuning as alternatives to the stoichiometric oxide systems [8]-[9]. A range of reported energy densities, benefits, and drawbacks can be found in Appendix A. 1.3.2. The Manganese Redox System and Magnesium-Manganese Oxide As mentioned, the manganese oxide compounds are interesting TCES candidates for their abundance and low toxicity to humans. The system is multi-valent, exhibiting several phase 8 transformations through reduction with increasing temperature. The sequence of reactions is given by Eq. 1.4, with temperature increasing from left to right. MnO2 → Mn2 O3 →Mn3 O4 →MnO (1.4). Of the three, the first reaction (MnO2/Mn2O3) is non-reversible and thus irrelevant for TCES [8]. The remaining two reactions are given by Eq. 1.5 and 1.6, along with their theoretical reaction enthalpies [8]. 6 Mn2 O3 ↔ 4 Mn3 O4 + O2 ΔHRxn = 202 kJ kg-1 (1.5). 1 Mn3 O4 ↔ 3MnO + 2 O2 ΔHRxn = 897 kJ kg-1 (1.6). Of these, the literature on the Mn2O3/Mn3O4 pair is far more expansive. The consensus is that this particular reaction is more conducive to TCES because it occurs below 1000 °C and has demonstrated reversibility, while the Mn3O4/MnO reaction occurs at very high temperatures (1500 °C) and has not yet demonstrated stable reversibility in the pure form. The literature on the Mn2O3/Mn3O4 reaction demonstrates poor oxidation kinetics, a susceptibility to sintering with accompanying reactive stability loss over many cycles, and a roughly 200 °C thermal hysteresis between the reduction temperature and the onset of oxidation [8]-[28]. Subsequent studies have focused on improving the energy density, reaction kinetics and cyclical stability within the 600- 1000 °C range by mixing Mn2O3 with other compounds (cobalt, iron, copper, lithium, and chromium), with the addition of iron proving to be the most successful. A collection of thermogravimetric, differential scanning calorimetry, and in-situ X-ray diffraction measurements on the Mn-Fe-O redox system over a range of iron concentrations says that compositions with around 20 mol% Fe have the best overall performance, reporting theoretical and experimental energy densities of 267 and 171-188 kJ kg-1, respectively, enhanced oxidation kinetics, a 9 reduction temperature of 995 °C and reactive stability over 75 cycles, with no noticeable effects from different synthesis methods [22]-[28]. Perhaps the most interesting of these when considering relevance to this work is the work by Wokon et al [27],[28]. An initial thermodynamic study for a material of Fe/Mn molar ratio of 1/3 used thermogravimetric analysis to report an energy density of 271 kJ kg-1 and demonstrated reactive stability over 100 cycles [27]. In a follow-up study, the group synthesized a 500 g sample of Fe-Mn-O of the same composition with the solid-state milling method from pure Fe2O3 and Mn3O4 powders, calcining in air at 1050 °C and sieving such that the final sample had a particle size range of 1-3mm. Lab-scale tests were performed with the particles in a packed bed tube reactor with an N2 and O2 dry-air-simulant mixture utilized as both heat transfer fluid and reactant. The system was cycled between 940-1040 °C and the oxygen content of the gas flow monitored so that the large-scale reduction extent could be compared with the small-scale from the thermogravimetric analysis. The group successfully demonstrated large-scale reactive stability over 17 cycles, achieving full reduction conversion after approx. 130 minutes of direct air heating [28]. It thereby seems reasonable to conclude that the Mn2O3/Mn3O4 couple is a viable energy storage material when mixed with iron on the order of 20-25 mol% Fe and is perhaps the most promising candidate material based on the Mn2O3/Mn3O4 redox pair. A much larger theoretical reaction enthalpy than that of the Mn2O3/Mn3O4 couple can be used if the Mn3O4/MnO redox pair can be made stable. Recall that the high reduction temperature of 1500 °C has been an obstacle to utilizing this reaction for TCES because of unproven reversibility and likely sintering via proximity of the reduction temperature to the melting temperature of Mn3O4 (1567 °C). To counteract this, Randhir et. al. [29] have demonstrated that the addition of magnesium oxide (MgO) to manganese oxide introduces 10 cyclical stability at temperatures between 1000-1500 °C because of the high melting temperature of MgO (2852 °C), reducing sintering and enhancing the structural stability. The detailed chemistry and stoichiometry are outlined in [29], with the simplified reaction given as y1 - y2 Mgx MnO1+x+y1 ↔ Mgx MnO1+x+y2 + 2 O2 (1.7) where x is the molar ratio of Mg to Mn and y1 and y2 are the excess oxygen present in the material at the oxidized and reduced states, respectively. Approx. 10 g samples of material with Mn/Mg molar ratios of 2/3, 1/1, and 2/1 are cycled in a lab-scale tube furnace reactor between 1000-1500 °C for 10 cycles under a constant air- simulant mixture of 80% argon and 20% oxygen. The oxygen concentration of the outlet gas is monitored using mass spectroscopy. The results show a cyclically stable reactivity (where “reactivity” refers to the total volume of oxygen released and absorbed per unit mass of the fully oxidized solid material over one full redox cycle) over 10 cycles, thus demonstrating short term reactive stability of the material. Cycling data for the three samples are shown in Figure 1.3. The cyclically stable reaction thus makes this material a promising candidate for a TCES reactor. 11 Figure 1.3. Reactive Cycling Stability. Reactive cycling stability of magnesium-manganese oxide of Mn/Mg molar ratios of (a) 2/3, (b) 1/1, and (c) 2/1 over 10 redox cycles [29]. Cycling is performed non-isothermally between an oxidation equilibrium state of (TOX = 1000 °C, PO2,OX = 0.2 atm) and a reduction equilibrium state of (TRED = 1500 °C, PO2,RED = 0.2 atm). Partial pressure of oxygen is held constant throughout the cycling. Data are presented as the volumetric release/absorption rate of oxygen gas per unit mass of fully oxidized solid material. (a) (b) (c) 12 1.3.3. Thermophysical Properties of Materials for Reactor Design For any novel TCES system, there are a few questions that must be answered: 1. Can the system reliably and repeatedly store and discharge energy, i.e., is it stable and reversible? 2. What is the maximum energy density of the system? 3. At what rates is energy discharged from the system? If a TCES system is to be designed around a packed bed of magnesium-manganese oxide, these questions need to be answered. Since reactive stability is already demonstrated in Figure 1.3, it is critical that the energy density (containing both sensible/latent and chemical reaction energy components) and the rate of the oxidation reaction are quantified for design of a magnesium-manganese oxide TCES reactor to proceed. An ideal redox metal oxide TCES reactor consists of a single vessel filled with a packed bed of reactive particles where air flows through the bed as both the reactant and the heat transfer fluid. In this system, the energy stored in the bed is directly transferred to the air as heat through two mechanisms: (1) particle-to-gas radiation and (2) interstitial convection. The overall heat transfer rate is a critical design point. This is highlighted by the oxidation step of the lab-scale study by Wokon et al. [28] of the reactive packed bed of the Fe/Mn of molar ratio 1/3 oxide material. In this step, the bed was slowly discharged from 1040 °C to 940 °C using the air- simulant gas flow. Temperature of the gas was precisely controlled at the inlet such that it was decreasing at a rate of 5 K/min for the duration of the discharge. From temperature data recorded at a few locations axially along the bed, they showed that a temperature front proceeded along the bed where the initial heat transfer was dominated by the removal of sensible heat from the system, followed by the onset of oxidation as the temperature locally was lowered [28]. This 13 demonstrates that the rate of energy removal from the reactor depends on the coupled process of interstitial convection, gas-particle radiation, and chemical oxidation. Depending on the oxidation kinetics of the material in question, the interstitial convective heat transfer could play a critical role in discharge rates and thus in reactor design and system control. These processes need to be independently quantified and then coupled with the oxidation kinetics in a system- level model for thorough and rigorous reactor design. It is important to note that detailed modeling of radiation in porous materials is very complex and is beyond the scope of this work. In a simplified approach, the radiative effects within a packed bed are accounted for by high-temperature effective thermal conductivity measurements taken for a stagnant packed bed of particles. This work was done for different beds of magnesium-manganese oxide particles in a previous study [30] with the results incorporated into this study to account for both particle-particle conduction within the bed and particle-gas radiative heat transfer. 1.4. Objectives and Organization The goal of this dissertation is to develop a simple coupled model of the thermodynamics, oxidation kinetics, and interstitial convective heat transfer for use in the design of a TCES reactor based on a packed bed of magnesium-manganese oxide particles. The model will serve as a tool for quantifying the energy discharge performance of a given packed bed. In particular, this model will determine the conditions under which the interstitial convective heat transfer plays a significant role in the rate of energy removal when coupled with the oxidation reaction (if any). To do so, this work undertakes the following tasks: 1. Measure the total energy density of magnesium-manganese oxide. Specify the distribution between sensible/latent energy and chemical reaction energy storage. 14 2. Identify the optimal operating conditions for a magnesium-manganese oxide TCES system for maximizing the total volumetric energy density (oxidation/reduction thermodynamic equilibrium states and specific material composition). 3. Quantify the reaction energy release rate with a bulk oxidation kinetic model for a packed bed of magnesium-manganese oxide particles. 4. Quantify the interstitial convective heat transfer for packed beds of magnesium- manganese oxide particles. 5. Model the energy discharge step of some different cases (inlet Reynolds number and packed bed aspect ratio) for a reactive magnesium-manganese oxide packed bed to determine which situations will experience a significant interstitial convection effect and which will not. A review of thermochemical energy storage, candidate materials, specific redox metal oxide materials, a discussion on the viability of magnesium-manganese oxide, and the work needed to model the coupled oxidation and interstitial heat transfer effects within a reactive packed bed of magnesium-manganese oxide particles has been discussed. The remainder of this dissertation is organized as follows: • Chapter 2 focuses on measuring the thermodynamic properties of magnesium-manganese oxide (total enthalpies, standard enthalpies of formation) needed to compute energy densities for the redox cycle operating between the equilibrium states of (TRED = 1500 °C, PO2,RED = 0.2 atm) and (TOX = 1000 °C, PO2,OX = 0.2 atm). The optimal thermodynamic equilibrium states for maximizing volumetric energy density are considered, and bulk oxidation rates are measured. 15 • Chapter 3 provides an overview of the aspects of porous media heat and mass transfer theory relevant to the problem of single-phase, non-thermal equilibrium flow through an isotropic, homogeneous porous medium and gives a background for the modeling approach used in Chapters 4 and 5. • Chapter 4 reviews the literature of experimental measurements of the interstitial heat transfer coefficient for packed beds of particles, details the experimental methods used in this work, and provides the results for the interstitial convective heat transfer in packed beds of magnesium-manganese oxide. • Chapter 5 couples the thermodynamic and kinetic measurements from Chapter 2 with the interstitial convective heat transfer measurements from Chapter 4 to build a simple coupled reactive packed bed model. A few cases are simulated for different inlet Reynolds numbers and bed length-to-radius (L/R) ratios to see when the interstitial convective heat transfer has a significant effect on the outlet behavior of a reactive packed bed of magnesium-manganese oxide. 16 CHAPTER 2: THERMOCHEMICAL PROPERTIES OF MAGNESIUM-MANGANESE OXIDE1 2.1. Introduction The first piece of thermodynamic information needed to model a novel TCES system is the total energy density of the storage material. For a redox metal oxide, this depends on the thermophysical properties (total enthalpies), the energy stored within the crystal lattice as the material is reduced (dependent on the standard enthalpies of formation of the oxidized and reduced phases), and the thermodynamic equilibrium states between which the system operates (oxidation/reduction temperatures, the partial pressure of oxygen during oxidation/reduction, and the specific molar composition of the material). Quantifying the total energy density is thus achieved by measuring the total enthalpy and standard enthalpies of formation of the oxidized and reduced phases at the thermodynamic equilibrium states between which the TCES system operates. The second piece of thermodynamic information needed is the optimal set of thermodynamic equilibrium states for bounding the thermochemical cycle to maximize the total energy density. For a redox metal oxide, this is achieved by measuring energy density and chemical reactivity of the material for different sets of thermodynamic equilibrium states. In the case of a fixed temperature cycle, this condenses to (1) different partial pressures of oxygen and (2) different chemical compositions of the metal oxide. In the case of a packed particle bed within a single reactor vessel, the limiting factor for how much energy can be stored is the 1 This chapter contains work previously published in [29],[31]-[33] in an exact, modified, or summarized form with permission from the journal editors. All work contained in this chapter was completed as a collaborative effort with Dr. Kelvin Randhir. 17 internal volume of the reactor vessel. Thus, the optimal set of equilibrium states will be those that maximize energy density on a per unit volume basis. Finally, the third piece of thermodynamic information needed is the rate at which the materials release the reaction energy. For a redox metal oxide, this involves quantifying the oxidation kinetics between the chosen equilibrium bounds of the thermochemical cycle. This chapter details the work done to measure the important thermodynamics of magnesium-manganese oxide operating under a thermochemical cycle between the oxidation and reduction equilibrium temperatures of 1000 and 1500 °C, respectively. Theoretical reactivity and energy density predictions for a few compositions are made using CALPHAD modeling. Total enthalpies and standard enthalpies of formation for the oxidized and reduced states are measured using a combination of drop and acid solution calorimetry. Total energy densities between 1000- 1500 ºC for the examined compositions are computed from these measurements. The extent to which energy density can be increased by lowering the partial pressure of oxygen during reduction is examined. Further compositional variation is studied for different molar ratios of Mn/Mg. From this work, a specific set of thermodynamic equilibrium states are suggested as operational bounds for maximizing the volumetric energy density of the system. Finally, the rate of reaction energy release is studied for a few Mn/Mg molar ratios by performing isothermal bulk oxidation kinetic measurements at temperatures between 1000-1500 ºC, from which a bulk oxidation kinetic model is derived. 2.2. Predictions from CALPHAD Models The Calculation of Phase Diagrams (CALPHAD) method is a computational thermodynamic tool used to model the equilibrium phase diagram for a given chemical system across a range of temperatures, pressures, and compositions [34]. Constituents can consist of 18 elements, vacancies, molecules, or ions in a disordered or ordered structure. Each constituent has its own temperature-dependent Gibbs free energy expression consisting of the component free energies: the reference state, the contribution from the atomic/molecular structural configuration (crystal lattice for ordered structures), the contributions due to physical phenomena, and any deviations from the previous three. Compositional variation is described using the Compound Energy Formalism [35] to model the interactions between all the constituents on different sublattices. Effects of gas pressure on the equilibrium are also accounted for in any system where gaseous components can exist. For a given set of thermodynamic state variables, a CALPHAD model repeatedly sums the different possible component Gibbs free energies for all potential phases until the minimum free energy is found. The composition associated with the minimum Gibbs free energy is thus the equilibrium state of the system for the given state variables. More discussion of CALPHAD modeling and the Compound Energy Formalism can be found in [34],[35] . CALPHAD modeling is a useful tool for predicting the reactivity and energy density of a material between two equilibrium states of interest. At the time of writing, two CALPHAD models exist in the literature for the Mg-Mn-O system. Panda et al [36] collected the available literature data for components of the Mg-Mn-O system and calculated an optimized set of model parameters using FactSage thermochemical software. They were able to reproduce the literature experimental data within their given error limits for the MgO-MnO-Mn2O3 system, along with categorizing the MgMn2O4 spinel structure. Dilner et al [37] performed a similar study using a different collection of literature data. They also performed specific heat measurements, differential thermal analysis (DTA), and X-ray diffraction (XRD) scans of three Mg-Mn-O 19 samples with different compositions to add data for the MgMn2O4 spinel and Mg6MnO8 phases and then calculated their own set of CALPHAD parameters. MATLAB programs were written using both models and used to predict the reactivity and the corresponding energy density for the 2/3, 1/1, and 2/1 Mn/Mg materials between 1000- 1500 °C for a few oxygen partial pressures of reduction. Reactivity data (i.e., the volume of oxygen released/absorbed per unit mass of the fully oxidized solid material over one full redox cycle) during stable cycling between 1000-1500 °C had already been collected for the 2/3, 1/1, and 2/1 Mn/Mg materials [29]. When the predictions from each model were compared with the cycling data, the Panda et al [36] model proved to be more accurate at high temperatures than the Dilner et al [37] model and was thus selected for energy density predictions. Figure 2.1 shows the predicted reactivity for 2/3, 1/1, and 2/1 Mn/Mg under oxygen partial pressures of reduction of 0.2 and 0.01 atm, along with the reactivity measured from the cycling data in Figure 1.3. Figure 2.2 predicts the energy densities of the same samples under oxygen partial pressures of reduction of 0.2 and 0.01 atm. Figures 2.1 and 2.2 are relative to a baseline oxidation equilibrium state of (TOX = 1000°C, PO2,OX = 0.2 atm). A detailed explanation of the error bars for the reactivity measurements can be found in Appendix B. 20 Figure 2.1. Reactivity Predictions. Reactivity of magnesium-manganese oxide (volume of oxygen released/absorbed per unit mass of fully oxidized solid material) with Mn/Mg molar ratios of 2/3, 1/1, and 2/1 assuming a reference oxidation equilibrium state of (TOX = 1000°C, PO2,OX = 0.2 atm) [29]. The CALPHAD model of Panda et al [36] predicts the reactivities at reduction equilibriums of (TRED = 1500 °C, PO2,RED = 0.2 atm) and (TRED = 1500 °C, PO2,RED = 0.01 atm). Measured reactivities from stable cycling between the reference oxidation equilibrium state of (TOX = 1000°C, PO2,OX = 0.2 atm) and the reduction equilibrium state of (TRED = 1500 °C, PO2,RED = 0.2 atm) are also included for comparison. CALPHAD predictions are in reasonable agreement for the 2/3 and 2/1 variants, but overpredict for the 1/1 variant. The model also predicts significant gains in reactivity by lowering PO2,RED, particularly as the Mn-content increases. 60 CALPHAD, 0.2 atm 50 Experimental, 0.2 atm CALPHAD, 0.01 atm 40 VO2 (cm 3 g -1) 30 20 10 0 2/3 1/1 2/1 Mn/Mg 21 Figure 2.2. Energy Density Predictions. Energy densities of magnesium-manganese oxide predicted from the CALPHAD model for molar ratios of Mn/Mg of 2/3, 1/1, and 2/1 for cycling between reduction equilibrium states of (TRED = 1500 °C, PO2,RED = 0.2 atm) and (TRED = 1500 °C, PO2,RED = 0.01 atm) are shown, assuming a reference oxidation equilibrium state of (TOX = 1000 °C, PO2,OX = 0.2 atm) [29],[36]. The higher reactivities predicted by the model in Figure 2.1 for PO2,RED = 0.01 atm correspond to higher predicted energy densities. 1400 0.2 atm 1200 0.01 atm 1000 ΔHS (kJ kg -1) 800 600 400 200 0 2/3 1/1 2/1 Mn/Mg The figures show that for the 1000-1500 °C thermochemical cycle, the Panda et al [36] model predicts energy densities greater than the theoretical energy storage density of the pure Mn2O3/MnO redox couple of 897 kJ kg-1 [8],[29]. However, the measured reactivity data are somewhat inconsistent from the CALPHAD model predictions, suggesting that experimental validation of the energy density is necessary. Thus, the following sections detail the measurement of the relevant thermodynamic properties (total enthalpies, standard enthalpies of 22 formation of the oxidized and reduced phases) and subsequent calculation of the energy storage densities of magnesium-manganese oxide. 2.3. Calorimetric Method for Measuring Energy Density 2.3.1. Introduction to Calorimetry with High-Temperature Metal Oxides The only way to measure thermodynamic properties of materials is by calorimetry, the science of measuring heat. Calorimetry is defined by the American Heritage Dictionary of the English Language as “the measurement of the amount of heat released or absorbed in a chemical reaction, change of state, or formation of a solution.” By using a calorimetric approach, the change in the state functions of a material can be determined between two given thermodynamic equilibrium states and is thus the correct experimental approach for measuring the thermodynamic properties of magnesium-manganese oxide. A wide variety of calorimetric techniques exist in the literature. The most common technique for measuring thermodynamic properties of high-temperature thermochemical materials is differential scanning calorimetry (DSC). DSC operates by gradually increasing the temperature of a sample and a reference material with a well-known heat capacity while maintaining both at the same temperature relative to each other. The difference in heat flows to the sample and reference is recorded and used to quantify the amount of energy absorbed/released during a chemical reaction or phase transition. Specific heat can also be measured using this approach. DSC is widely used to quantify the energy densities of the redox metal oxides reviewed in Chapter 1 [8]- [28] . DSC is a useful tool, but has some notable drawbacks when applied to high-temperature [38] redox metal oxides. Jacob et al examined the accuracy and reproducibility of specific heat measurements using DSC at high temperatures (> 700 °C), finding that results at temperatures ≥ 23 800 °C were highly variable with variations in sample mass, mass of the calibration standard, heating rate, and flow rate of the purge gas. Saeed et al [39] studied the uncertainties associated with the enthalpy of melting and transition temperature measurements of phase change materials taken using DSC, finding that high heating rates and/or large sample masses resulted in shifts towards results of higher magnitude. In their work, the low thermal conductivities and high energy densities of phase change materials created a noticeable temperature gradient within the sample, causing the DSC control thermocouple to measure a higher temperature than that of the actual average sample temperature, creating a measurement bias. A smaller sample mass and slower heating rate reduced the thermal gradient bias, but widened the range of uncertainty because of a low observable heat flux [39]. These two studies hold relevance for redox metal oxides at high temperatures (≥ 800 °C). The first implies that the low accuracy and variable reproducibility above 800 °C will increase the uncertainty in any enthalpy measurements collected using DSC at the temperatures of interest for redox metal oxide TCES. In addition, redox metal oxides have high energy storage densities and low thermal conductivities (when subjected to high temperatures). Therefore, the temperature bias caused by the required fast heating rate and inherent temperature gradient could make it difficult to accurately measure the reaction enthalpies of redox metal oxides at high temperature using DSC. Previous attempts at magnesium-manganese oxide energy density measurements using DSC were not satisfactorily accurate due to these measurement issues and thus directed the development of a different calorimetric approach to improve the accuracy. 2.3.2. Calorimetric Approach To this end, a two-part approach was developed using two different calorimetry techniques: drop calorimetry and acid solution calorimetry. Drop calorimetry measures the total enthalpy 24 change of a material between two temperatures by bringing a sample to thermal equilibrium at some higher temperature, then rapidly cooling it by dropping it into a well-insulated device at some lower reference temperature (generally room temperature). The change in temperature of the calorimeter is recorded until thermal equilibrium is reached. From this, the total enthalpy of the sample at the drop temperature can be calculated relative to the reference temperature. Measurements at a series of drop temperatures facilitate the measurement of the specific heat of the material. An effective drop calorimeter will prevent any reaction or phase change of the sample such that the chemical state is known, be capable of accurately recording the heat released by the sample, and have a quantifiable heat loss such that it can be accounted for in the measurement [40]. Acid solution calorimetry is a type of calorimetry in which a sample is dissolved in an acidic solution. The dissolution breaks the chemical bonds and releases the bond-energy into the solution, creating a measurable temperature rise. When starting with a room temperature solution, this approach measures the standard enthalpy of formation of phases present in the sample when dissolved: by recording the change in solution temperature from the initial to the final state, one can use a series of energy equations derived using Hess’s law from the collection of chemical reactions that construct the chemical system from the base elements to quantify the bond-energy of the phases present in the sample. By combining these two approaches, one can measure the enthalpies necessary for computing the energy density of a material between the oxidized and reduced equilibrium states. Both techniques have been successfully utilized in the literature. Drop calorimetry has frequently been used for measuring total enthalpy of metal oxides at high temperatures. For example, Fredrickson et al [41] measured the enthalpy of uranium dioxide up to 1227 °C. Popa et al [42] measured the thermodynamic properties of monazite from 177-1227 °C. Takahashi et al [43] 25 measured the enthalpy and specific heat of gadolinia-doped uranium dioxide between 127-727 °C. Ritchet et al [44] measured the thermodynamic properties of several silicon dioxide phases between 727-1527 °C. Acid solution calorimetry has also been shown to be a successful technique. Brink and Holley [45] used this technique to calculate the standard enthalpy of formation of strontium (II) oxide by dissolving it and pure strontium metal in aqueous 1 M hydrochloric acid solution (HCl). [46] Matskevich et al successfully calculated the standard enthalpy of formation of SrCeO3 perovskite doped with lutetium (III) oxide by dissolving samples in a solution of 1 M HCl. [47] Tsvetkov et al calculated the standard enthalpy of formation of double perovskites GdBaCo2−xMxO6-δ (M= Fe, Mn; x = 0, 0.2) by dissolving constituent metal-oxides in 200 mL of 4 M HCl solution. Figure 2.3 illustrates how these two techniques are applied to determine the energy density of magnesium-manganese oxide. The two techniques break the energy density down into two separate measurements of the sensible/latent and chemical reaction energies. During the drop step, cooling the material under an inert gas allows for the latent and sensible heat to be released while preventing heat release due to oxidation. By doing so, the enthalpy at temperatures of interest can be measured for both the oxidized and reduced states, allowing the isolation of the energy stored as sensible and latent heat as the temperature is increased from the oxidation to the reduction equilibrium temperature. The chemical reaction energy can then be calculated as the difference between the standard enthalpies of formation of the reduced and oxidized states from the acid solution step. The total energy stored for a specified temperature cycle is then computed by, ΔHS = ΔHDrop, RED – ( ΔHDrop, OX – ΔHRxn) (2.1) 26 where ΔHDrop, RED is the enthalpy difference between 1500 and 25 °C for reduced materials, ΔHDrop, OX is the enthalpy difference between 1000 and 25 °C for oxidized materials, ΔHRxn is the chemical energy stored during thermal reduction, and ΔHS is the total energy density. Figure 2.3. Calorimetric Method Diagram. Hypothetical enthalpy curves for the oxidized and reduced materials. ΔHDrop, RED and ΔHDrop, OX are the total enthalpies of the reduced and oxidized phases, respectively, measured from drop calorimetry. The ΔHf terms are the standard enthalpies of formation of the oxidized and reduced phases and are measured from acid solution calorimetry. y is the excess oxygen beyond the stoichiometric present in the simplified chemical formula for magnesium-manganese oxide. ΔHS is the total energy density for the thermochemical cycle between 1000-1500 °C and is calculated from Eq. 2.1. Both curves are for constant PO2. To obtain accurate acid solution calorimetry measurements, the material must dissolve quickly within the acid. When the material takes a long time (order of hours) to dissolve, the concurrent heat loss from the system reduces the accuracy of the calorimetry data. Many metal oxides, in particular, dissolve very slowly in strong acid solutions. However, by adding the right 27 catalyst to the acid, the rate of dissolution of the material can be enhanced enough that accurate measurements can be taken. For the purposes of this work, adding tin (II) chloride (SnCl2) to aqueous hydrochloric acid solution enhanced the rate of dissolution of magnesium-manganese oxide sufficiently for acid solution calorimetry. The experimental setups, procedure for the two- part method, calorimeter calibration/validation, and justification for the addition of SnCl2 are described in the following sections. 2.4. Experimental Setup and Procedure 2.4.1. The Drop Calorimeter A schematic diagram of the drop calorimeter used in this work is shown in Figure 2.4. The system consists of a 44.45 mm OD stainless steel tube with height 152.4 mm that is closed at one end and placed inside a plastic Styrofoam-lined thermos. 316 stainless steel is used to prevent reaction of the container with dropped samples. Copper fins are attached to the exterior of the tube to enhance the heat transfer to the surrounding fluid. The thermos is filled with 780 mL of distilled water (specific heat of approx. 4.2 kJ kg-1 K-1 at room temperature and boiling point of approximately 100 °C) to serve as the calorimetry medium. The temperature rise of the water is recorded using K-type thermocouples. Boiling of the water is of no concern, as only small increases from room temperature occur. Argon is routed to flow through the tube apparatus to purge air and maintain a minimal oxygen partial pressure. A Styrofoam lid is then placed over the thermos to fully insulate the system, with a hole cut in the center through which the sample is dropped. A Styrofoam plug seals the hole during the quenching process. During pre-heating, the sample is suspended inside a 50.8 mm OD alumina tube mounted within a Sentrotech high-temperature furnace capable of reaching 1700 ˚C. A B-type ceramic- encased thermocouple is placed inside the tube within the heated zone of the furnace so that the 28 exact temperature of the sample can be monitored and precisely controlled. High temperature insulation from Zircar Zirconia is used as a heat barrier inside the end of the tube. Prior to the drop experiments, the insulation is heat-treated for 4 hours at 600 °C to burn off the organic binding that can cause the insulation to stick to the alumina tube and thereby avoid any difficulty when trying to remove it. Figure 2.4. Drop Calorimeter. Schematic diagram of the drop calorimeter used to measure the total enthalpies of magnesium-manganese oxide at high temperatures [31]. Total enthalpies are relative to room temperature. Argon gas is routed through the stainless steel tube to impose an inert environment during drop experiments to avoid any oxidation of the sample. 2.4.2 The Acid Solution Calorimeter A schematic diagram of the acid solution calorimeter is shown in Figure 2.5. The system consists of approximately 200 g of 5.45 M hydrochloric acid solution contained within a Dewar flask. This concentration was chosen because its high boiling point of approximately 108 °C (relative to stronger concentrations) would prevent boiling due to potentially high instantaneous heat release during dissolution and avoid error from latent heat loss. A foam insulating cover is 29 tightly fixed inside the top of the Dewar flask to minimize heat loss through the top. PFA (perfluoro alkoxy) coated K-type thermocouples are used to measure temperature. PFA coating is necessary because it prevents reaction of the thermocouple with the acid. A magnetic stirrer is placed inside the flask to agitate the solution and increase the dissolution rate. Stirring the pure solution with the magnetic stirrer for an extended length of time did not create any measurable temperature rise within the solution and was thus deemed reliable to use for experiments. Figure 2.5. Acid Solution Calorimeter. Schematic diagram of the acid solution calorimeter [31] . 2.4.3. The Calorimetric Procedure To begin the measurement method, samples are prepared by sintering a mixture of the raw component oxides, mixed to the desired molar ratios, at 1200 °C under 0.1 SLPM of hydrogen and 0.4 SLPM of argon for approximately 12 hours to form solid structures of the pure component oxides. Approximately 2-4 g of the prepared sample is placed within a small alumina ceramic crucible. The initial mass is recorded, as it is needed for the acid solution calorimetry calculations. 30 This sample preparation step is necessary because the value of y in the simplified chemical formula for magnesium-manganese oxide (Eq. 1.7) is defined as the excess oxygen present in the crystal lattice beyond that of the pure component oxides (MnO and MgO). The y values are later calculated as the difference between this initial sample mass and the mass of the quenched sample for both reduced and oxidized materials. The sample is suspended inside the tube furnace, heated to the reduction temperature, and allowed to dwell for approx. 4 hours. The sample must be kept at a constant oxygen partial pressure by means of constant gas flow rates through the tube. For this work, flow rates of 0.8 and 0.2 SLPM of argon and oxygen, respectively, are chosen to maintain an oxygen partial pressure of approximately 0.2 atm. If studying reduced materials, the sample is ready to be quenched after the prescribed dwell time. If studying oxidized materials, the sample is then lowered to the oxidation temperature of interest and allowed to dwell for an additional 4 hours to ensure chemical equilibrium is achieved. After the dwell, the drop calorimeter is flooded with a high flow rate (approximately 6 SLPM) of argon. This is crucial for preventing any oxidation of the sample to ensure that the phase structure that exists at the dwell temperature is maintained. If not, the material [36], [37] would partially oxidize during the quench to form phases that exist at low temperatures , altering the enthalpy measurement. The insulation from the bottom end of the tube inside the furnace is removed to clear a path for the sample to drop. The calorimeter is placed as close as possible to the mouth of the tube, and the sample is dislodged and rapidly quenched inside the calorimeter. The lid of the calorimeter is immediately plugged to isolate the system and prevent heat loss. The heat of the sample is dissipated to the water surrounding the stainless steel tube inside the device. The temperature of the water is recorded, and the device remains quiescent until the water reaches a constant temperature. Following the temperature measurement, the sample is 31 removed and the final mass is measured. This mass is the mass of the sample after undergoing transformation from the pure oxides to the phases that exist at the drop temperature and is defined as the final mass. The sample is then ground in a mortar and pestle to a fine powder under acetone to prevent oxidation during grinding and thoroughly dried at a low temperature inside a box furnace (around 50 °C) or by forced convection with a fan. After drying, the mass of acid is measured in the Dewar flask and approximately 5.5 g of SnCl2 is added. The acid is kept at room temperature and is stirred inside the flask for 5-10 minutes to allow it to reach thermal equilibrium with the walls of the flask and to fully dissolve the SnCl2. The thermocouple is inserted into the solution through the lid, securely fitted at the top to minimize heat loss. The initial temperature of the solution is recorded. A small sample of finely-ground material is then dissolved in the solution, the lid securely reaffixed after introduction of the sample, and the temperature rise is recorded. Upon total dissolution, the temperature will reach a constant value and is recorded as the final temperature. The contents of the flask are emptied and the flask is cleaned. The experiment is performed for varying sample masses ranging between 0.1 and 2 g, for a total of 3-4 trials. Specific masses were chosen from within this range based on the availability of material and on the speed of dissolution of the oxides in the prescribed volume of acid, with larger masses chosen for oxides with faster dissolution rates. The resulting temperature rises are plotted versus sample mass to generate a linear curve fit through the origin. The results of the experiments are used to calculate the total energy density of the material at the drop temperature relative to 25 °C, as well as to demarcate the portion of the total stored as sensible/latent and that stored as chemical reaction energy. This procedure is repeated for other thermodynamic equilibrium states of interest so that the energy density for a given redox cycle can be computed. 32 2.5. Calorimeter Calibration and Validation 2.5.1. The Drop Calorimeter Calibration of the drop calorimeter was done with well-studied, non-reactive materials with known heat capacities. Aluminum (III) oxide (α-Al2O3) of 99.9% purity purchased from Sigma- Aldrich was chosen as the calibration standard. Thermodynamic data were obtained for the enthalpy change of Al2O3 from 25 °C to the temperatures of interest [48]. Because calibration was performed with non-reactive materials, high-temperature dwell times from the general methodology were reduced from 4-8 hours to 1 hour as only thermal equilibrium needed to be achieved, after which samples were then rapidly quenched in the calorimeter. This was repeated for several temperatures ranging from 1000-1500 °C. The known changes in enthalpy were plotted versus the recorded temperature rises and a linear fit through the origin was generated, taking the calibration constant as the slope of the fit, given by ΔHAl2 O3 C1 = ΔTH2 O (2.2). Verification of the constant was achieved by performing the same experiments with samples of zirconium (IV) oxide (ZrO2) of similar mass and purity of 99.5%, obtained from Stanford Materials Corporation. 2.5.2. The Acid Solution Calorimeter 2.5.2.1. Addition of SnCl2 to the Calorimetry Process When first attempting to dissolve magnesium-manganese oxide in hydrochloric acid, the dissolution rate was extremely slow. This presented an obstacle for measuring the standard enthalpies of formation of magnesium-manganese oxide phases using this method. Subsequent trials of dissolving pure magnesium (II) oxide (MgO), manganese (II) oxide (MnO), and manganese (III) oxide (Mn2O3) in hydrochloric acid suggested that the slow dissolution was due 33 to the presence of Mn3+ ions in the magnesium-manganese oxide. Initially, this hypothesis was simply observational: pure MgO and MnO, both containing a metal ion with +2 charge, dissolved quite rapidly, but Mn2O3, containing Mn3+ ions, did not dissolve well. Review of the literature supported this hypothesis. Parida and Das [49] showed that adding SnCl2 to hydrochloric acid will enhance the rate of dissolution of iron (III) oxide (Fe2O3). To see if the thermodynamics supported the use of SnCl2 for enhancing the dissolution rate of Mn2O3, standard reduction potentials were obtained from Bard et al. [50] for the reduction of Fe3+, reduction of Mn3+, and oxidation of Sn2+, shown in Table 2.1. Table 2.1. Standard Reduction Potentials. Standard reduction potentials for relevant metal ions [50].The smaller magnitude of Eo for half-reaction 1 relative to half-reactions 2 and 3 suggests that the incorporation of Sn2+ ions into the solution will reduce the Me3+ ions to Me2+ ions as their reaction potentials indicate a higher tendency to be reduced than the Sn2+ ions. Half-Reaction # Half-Reaction Eo (V) 1 Sn4+ (aq) + 2e- → Sn2+ (aq) 0.15 2 Fe3+ (aq) + e- → Fe2+ (aq) 0.77 3 Mn3+ (aq) + e- → Mn2+ (aq) 1.5 The Gibbs free energies of reaction for the dissolution of Fe2O3 and Mn2O3 in hydrochloric acid [50] with and without SnCl2 were also calculated . Table 2.2 gives the dissolution reactions, and Figure 2.6 shows the Gibbs free energy of reaction for each reaction in Table 2.2. 34 Table 2.2. Free Energy of Dissolution of Me3+ Compounds. Dissolution reactions of Fe2O3 and Mn2O3 in hydrochloric acid both without and with the addition of SnCl2 to the solution [50] , with reactions labeled by number for reference. Reaction # Dissolution Reaction ΔGRxn 1 Fe2O3 (s) + 6H+ (aq) → 2Fe3+ (aq) + 3H2O (l) ΔG1 2 Fe2O3 (s) + 6H+ (aq) + Sn2+ (aq) → 2Fe2+ (aq) + Sn4+ (aq) + 3H2O (l) ΔG2 3 Mn2O3 (s) + 6H + (aq) + 2Cl- (aq) → 2Mn2+ (aq) + 3H2O (l) + Cl2 (g) ΔG3 4 Mn2O3 (s) + 6H+ (aq) + Sn2+ (aq) → 2Mn2+ (aq) +3H2O (l) + Sn4+ (aq) ΔG4 Figure 2.6. SnCl2 Dissolution Enhancement, Predictions. The calculated Gibbs free energies of reaction for the dissolution reactions in Table 2.2 are shown. The values for the dissolution reactions with SnCl2 in solution are significantly more negative than those for dissolution reactions without SnCl2 in solution. This implies that chemical equilibrium for dissolution of the Me3+ ions is far more favorable for solutions with SnCl2, suggesting that adding SnCl2 could enhance dissolution rates of the Me3+ ion compounds. 35 The thermodynamics suggest that SnCl2 will enhance the dissolution rate of Mn2O3 in hydrochloric acid. From Table 2.1, we see that Eo2 > Eo1, indicating that the presence of Sn2+ ions in the solution will reduce Fe3+ to Fe2+. Further, ΔG2 < ΔG1, indicating that the thermodynamics are more favorable for dissolution of Fe2O3 when SnCl2 is included. These results are somewhat obvious, as it is already known that SnCl2 enhances the dissolution rate of Fe2O3 in hydrochloric [49] acid . However, these results are important for drawing comparisons between the Fe2O3 data and that for Mn2O3. The same trends are present for Mn2O3: Eo3 > Eo1, indicating that the presence of Sn2+ ions in the solution will reduce Mn3+ to Mn2+, and ΔG4 < ΔG3, showing that the thermodynamics are more favorable for dissolution of Mn2O3 with SnCl2 included in the solution. Therefore, it is hypothesized that the addition of SnCl2 will enhance the dissolution rate of Mn2O3 in hydrochloric acid and thereby will also enhance dissolution of magnesium-manganese oxide. To verify this, tests were conducted where approximately 1 g each of MnO, MgO, Mn2O3, and 2/1 Mn/Mg were dissolved for precisely 10 min in the aqueous hydrochloric acid solution. Upon completion of the 10 min, portions of the solutions were separated into vials and pictures were taken to show the extent of dissolution. This procedure was repeated for solutions with approximately 5.5 g of SnCl2 dissolved prior to addition of the oxides. The images showing the extent of dissolution can be seen in Figure 2.7. 36 Figure 2.7. SnCl2 Dissolution Enhancement, Test. From left to right, samples of MnO, MgO, Mn2O3, and 2/1 Mn/Mg after 10 minutes of dissolution in 5.45 M hydrochloric acid solution (a) without and (b) with SnCl2 in the solution. Filtration of the samples yielded noticeable levels of particulate for the Mn2O3 and 2/1 Mn/Mg samples without SnCl2 in solution. Conversely, filtration of the solutions with SnCl2 yielded no particulate matter for all four samples [31]. (a) Solution without SnCl2 MnO MgO Mn2O3 2/1 Mn/Mg (b) Solution with SnCl2 MnO MgO Mn2O3 2/1 Mn/Mg 37 In Figure 2.7(a) the two vials containing Mn2O3 and 2/1 Mn/Mg are both opaque solutions. Subsequent filtering of the solutions showed that a large amount of particulate remained undissolved. Meanwhile, the amount of particulate caught by the filtering process for the solutions containing MnO and MgO was negligible. However, in Figure 2.7(b), the vials containing Mn2O3 and 2/1 Mn/Mg are perfectly clear solutions, and subsequent filtering left no observable particulate behind for any of the four solutions. It is clear from these results that the addition of SnCl2 will enhance the rate of dissolution of the Mn3+ ions present in Mn2O3 and in magnesium-manganese oxide. The next step is to verify that the addition of SnCl2 to the aqueous hydrochloric acid solution enables acquisition of accurate acid solution calorimetry data for these higher valence state transition metal oxides. 2.5.2.2. Calibration & Validation The previously described methodology is used for calibrating the acid solution calorimeter and validating the method. MgO of 97%+ purity from Acros Organics is used as the calibration standard. MnO of 99% purity from Strem Chemicals and zinc (II) oxide (ZnO) of 99% purity from Fischer Scientific are used to validate the calibration. Once validated, SnCl2 of purity ≥ 99% from Alfa Aesar is added to the system, and experiments are run with both Mn2O3 and Fe2O3 from Strem Chemicals of purity 99% and 99.8%, respectively, to verify both that SnCl2 enables accurate acid solution calorimetry data and that its addition does not interfere with the system in a significant way. If true, the pre-existing calibration constant from the +2 metal ion compounds will still give accurate results for the higher valence state metal oxides. All chemical equations used for acid solution calorimetry in this work are defined in Table 2.3, along with their corresponding enthalpies of reaction. 38 Table 2.3. Acid Solution Calorimetry Chemical Reactions. The theoretical chemical reactions used throughout this work and the corresponding enthalpies of reaction. This collection of equations is used according to Hess’s Law to calculate standard enthalpies of formation of compounds dissolved in the acid solution calorimeter. Reaction Net Ionic Chemical Reaction ∆HRxn (kJ mol-1) # 1 1 Mg (s) + 2 O2 (g) → MgO (s) ΔH1 1 2 H2 (g) + 2 O2 (g) → H2O (l) ΔH2 3 Mg (s) + 2H+ (aq) → Mg2+ (aq) + H2 (g) ΔH3 4 MgO (s) + 2H+ (aq) → Mg2+ (aq) + H2O (l) ΔH4 1 5 Mn (s) + 2 O2 (g) → MnO (s) ΔH5 6 Mn (s) + 2H+ (aq) → Mn2+ (aq) + H2 (g) ΔH6 7 MnO (s) + 2H+ (aq) → Mn2+ (aq) + H2O (l) ΔH7 1 8 Zn (s) + 2 O2 (g) → ZnO (s) ∆H8 9 Zn (s) + 2H+ (aq) → Zn2+ (aq) + H2 (g) ∆H9 10 ZnO (s) + 2H+ (aq) → Zn2+ (aq) + H2O (l) ∆H10 3 11 2Mn (s) + 2 O2 (g) → Mn2O3 (s) ∆H11 # 1 12 $ H2 (g) + 2 Cl2 (g) → HCl (g) ∆H12 13 Sn (s) + Cl2 (g) → SnCl2 (s) ∆H13 14 Sn (s) + 2Cl2 (g) → SnCl4 (l) ∆H14 15 HCl (g) → H+ (aq) + Cl- (aq) ∆H15 16 SnCl2 (s) → Sn2+ (aq) + 2Cl- (aq) ∆H16 17 SnCl4 (l) → Sn4+ (aq) + 4Cl- (aq) ΔH17 39 Table 2.3 (cont’d). Mn2O3 (s) + 6H+ (aq) + Sn2+ (aq) → 2Mn2+ (aq) +3H2O (l) 18 ∆H18 + Sn4+ (aq) 3 19 2Fe (s) + 2 O2 (g) → Fe2O3 (s) ∆H19 20 Fe (s) + 2H+ (aq) → Fe2+ (aq) + H2 (g) ∆H20 Fe2O3 (s) + 6H+ (aq) + Sn2+ (aq) → 2Fe2+ (aq) +3H2O (l) + 21 ∆H21 4+ Sn (aq) (1+x+y) 22 xMg + Mn + 2 O2 → MgxMnO1+x+y ∆H22 MgxMnO1+x+y + 2(1+x+y)H+ + ySn2+ → xMg2+ + Mn2+ + 23 ∆H23 ySn4+ + (1+x+y)H2O The equations used to calculate acid solution calorimetry results are derived according to Hess’s law from the chemical equations provided in Table 2.3. It is first necessary to define the expression for calculating the enthalpy of dissolution of a compound in the acid (the quantity being directly measured by the calorimeter). This is given by, C2 M ΔHDiss = - m ΔT (2.3). Here, ΔHDiss is the enthalpy of the dissolution reaction, C2 is the system-dependent calibration factor, M is the molecular weight of the sample, ∆T is the measured temperature rise of the solution, and m is the sample mass. The negative sign denotes the exothermic nature of the reaction. To calculate C2 using MgO, the theoretical reaction scheme for its dissolution in hydrochloric acid is written starting from the base elements with equations taken from Table 2.3. The system is balanced, simplified, and reduced to a single equation in terms of reaction enthalpies. Solving for the standard enthalpy of formation of MgO yields, 40 ΔH1 = ΔH2 + ΔH3 – ΔH4 (2.4). Here, ∆H1 is the standard enthalpy of formation of MgO, ∆H2 is the standard enthalpy of formation of liquid water, and ∆H3 and ∆H4 are the enthalpies of the dissolution reactions of pure magnesium metal (Mg) and MgO, respectively. The derivation of Eq. 2.4 from the reactions taken from Table 2.3 is detailed in Appendix B as an example of the procedure for deriving all subsequent acid solution calorimetry equations. Applying Eq. 2.3 to ΔH3 and ΔH4, substituting into Eq. 2.4 and solving for the calibration constant yields, ΔH1 - ΔH2 C2 = - ΔT ΔT (2.5). MMg & m ' - MMgO & ' m MgO Mg For this system, the constant is calculated to be 662 ± 10 J °C-1. The data for generating the linear fits used in calibration and throughout the rest of this work can be found in Appendix B. The constant is validated using both MnO and ZnO, following the same procedure as MgO. The standard enthalpy of formation for MnO is calculated by, ΔH5 = ΔH2 + ΔH6 – ΔH7 (2.6). ΔH5 is the standard enthalpy of formation of MnO, and ΔH6 and ΔH7 are the enthalpy changes of the dissolution reactions of Mn and MnO, respectively. Similarly, the standard enthalpy of formation of ZnO is given by, ΔH8 = ΔH2 + ΔH9 – ΔH10 (2.7). ΔH8 is the standard enthalpy of formation of ZnO, and ΔH9 and ΔH10 are the enthalpy changes of the dissolution reactions of pure zinc (Zn) metal and ZnO, respectively. The next phase required the addition of SnCl2 to the system to validate that it will enhance the dissolution rate of Mn2O3 and Fe2O3 such that the acid solution calorimetry data is accurate. The standard enthalpy of formation of Mn2O3 is given by, 41 ∆H11 = 3ΔH2 - ∆H18 + 2ΔH6 - 2∆H12 -2∆H15 + ∆H14 - ∆H13 + ∆H17 - ∆H16. (2.8). Here, ∆H11 is the standard enthalpy of formation of Mn2O3; ∆H12, ∆H13, and ∆H14 are the standard enthalpies of formation of gaseous hydrogen chloride (HCl), SnCl2, and tin (IV) chloride (SnCl4), respectively; ∆H15 is the enthalpy of the dissolution reaction of gaseous HCl in water to form hydrochloric acid; and ∆H16, ∆H17, and ∆H18 are the enthalpies of the dissolution reactions in aqueous hydrochloric acid of SnCl2, SnCl4, and Mn2O3, respectively. ∆H17 is measured as approximately -90.6 ± 1.4 kJ mol-1, while ∆H16 is neglected as no noticeable temperature change occurs upon dissolution of SnCl2 in the solution. ΔH15 was calculated for 5.45 M acid as -12.6 kJ mol-1 [51],[52]. Similarly, the standard enthalpy of formation of Fe2O3 is given by, ∆H19 = 3∆H2 - ∆H21 + 2∆H20 - ∆H13 + ∆H14 - ∆H16 + ∆H17 - 2∆H12 - 2∆H15 (2.9). Here, ∆H19 is the standard enthalpy of formation of Fe2O3, ∆H20 is the enthalpy of the dissolution reaction of pure iron powder (Fe), and ∆H21 is the enthalpy of the dissolution reaction of Fe2O3. All of the validation results are shown in Figure 2.8. 42 Figure 2.8. Acid Solution Calorimeter Validation. Standard enthalpies of formation for ZnO, MnO, Mn2O3, and Fe2O3. Values are shown from both the accepted literature [48] and measurements from the acid solution calorimeter. Error is within a maximum of 2.5% of the accepted values for all acid solution calorimetry measurements shown. It is clear that the calibration is quite accurate, as the standard enthalpies of formation measurements in Figure 2.8 all have error < 2.5% from the accepted values. The error analysis for the acid solution calorimetry measurements is found in Appendix B. Further, because the results obtained using SnCl2 for Mn2O3 and Fe2O3 have a similarly high accuracy, it is safe to conclude that adding SnCl2 enhances the dissolution rates of these higher valence state transition metal oxides enough that acid solution calorimetry is accurate. Therefore, the system is considered well-calibrated and suitable for measuring the standard enthalpy of formation of magnesium-manganese oxide. 43 In regards to the dissolution experiments involving the pure metals, there are a few important factors to note. For pure Mg and pure Mn, dissolution is done using sample masses ≤ 0.1 g for Mg and ≤ 0.5 g for Mn, with the metal contained inside a 1.7 mL Posi-Click polypropylene micro-centrifuge tube with a slit cut on each side to slow the exposure of the pure metal to the acid. This is necessary because the heat generated by exposing an entire sample to the solution instantaneously causes a rapid reaction, raising the temperature of the acid such that some boiling occurs and latent heat is removed from the system by escaping hydrogen gas, leading to measurement error. Further, dissolving samples larger than the stated upper mass limits also causes this problem. Therefore, samples of masses no greater than the given limits are used. Slow exposure is done such that the instantaneous temperature never rises above 50 °C, ensuring all heat generated is dissipated throughout the solution. A hole is cut in the cap of the vial for insertion of the thermocouple. The thermocouple is then used to gradually expose the sample to the acid, thereby preventing significant vapor formation. Additionally, for pure Fe, it is necessary to remove all magnetic components from the system and allow dissolution to proceed either naturally or by manually agitating the solution, as use of a magnetic stirrer causes some induction heating in the system and introduces error in the final temperature measurement. 2.6. Measurements for Magnesium-Manganese Oxide The calorimetric method is applied to magnesium-manganese oxide samples with Mn/Mg molar ratios of 2/3, 1/1, and 2/1. The total energy density for this temperature cycle is calculated using Eq. 2.1. Because the samples were contained within alumina crucibles during the drop calorimetry process, the enthalpy change of Al2O3 is calculated for each drop temperature using [48] thermodynamic data and is then subtracted from the total enthalpy change measured by the calorimeter to yield the drop calorimetry results for magnesium-manganese oxide. 44 To verify that the argon quench successfully prevents oxidation of the sample and preserves the compounds that are present at the drop temperature, XRD data were collected for samples rapidly quenched from drop temperatures of 1000 and 1500 °C, shown in Figure 2.9. Figure 2.9. XRD Data of Quenched Samples. XRD data for samples quenched from (a) the reduction equilibrium state of (TRED = 1500 °C, PO2,RED = 0.2 atm) and (b) the oxidation equilibrium state of (TOX = 1000 °C, PO2,OX = 0.2 atm) [29]. The phases present in the reduced (primarily monoxide) and oxidized (primarily spinel) states agree well with the phase diagram of the Mg-Mn-O system calculated by Panda et al [36] and presented in their work. These results suggest that the quenching of samples in a high flow of inert gas inside the drop calorimeter preserves the phases that exist at the drop temperature. (a) 45 Figure 2.9 (cont’d.). (b) The XRD data show that the quenching process preserves the compounds that form at high temperatures, as the dominant peaks are the monoxide phase for reduced materials and the spinel [29] phase for oxidized materials, in accordance with the theory detailed by Randhir et al and the [36] Mg-Mn-O system phase diagram of Panda et al . Therefore, we can conclude that the methodology presented in this work is capable of accurate energy density measurements for TCES redox materials. The total enthalpies of magnesium-manganese oxide measured from the drop calorimetry process are shown in Figure 2.10. ΔHOX represents the total enthalpy change between 25 and 1000 °C for materials of the oxidized chemical formula MgxMnO1+x+y1, as detailed in Figure 2.3. Similarly, ΔHRED represents the total enthalpy change between 25 and 1500 °C for materials of the reduced chemical formula MgxMnO1+x+y2. 46 Figure 2.10. Drop Calorimetry Results. Total enthalpies of the oxidized and reduced magnesium-manganese oxide phases measured from drop calorimetry, relative to a room- temperature equilibrium state of (T = 25 °C, PO2 = 0.2 atm), for Mn/Mg molar ratio variants of 2/3, 1/1, and 2/1. The reduced phase is quenched from an equilibrium state of (TRED = 1500 °C, PO2,RED = 0.2 atm). The oxidized phase is quenched from an equilibrium state of (TOX = 1000 °C, PO2,OX = 0.2 atm) [31]. 1800 Oxidized Phase 1600 Reduced Phase 1400 1200 ΔHDrop (kJ kg-1) 1000 800 600 400 200 0 2/3 1/1 2/1 Mn/Mg For acid solution calorimetry, the standard enthalpy of formation of magnesium- manganese oxide is derived similarly to previous compounds from the relevant chemical reactions in Table 2.3 as ∆H22 = ∆H6 + x∆H3 + (1+x+y)∆H2 + y(∆H14 - ∆H13 + ∆H17 - ∆H16 - 2∆H12 - 2∆H15) - ∆H23 (2.10). Here, ∆H22 is the standard enthalpy of formation of the material, ∆H23 is the enthalpy of the dissolution reaction of the material, x is the molar ratio of Mg/Mn, and y is the excess oxygen atoms present in the material, with the unit moles of O per mole of material. y is calculated as the 47 difference between the final mass of the quenched sample and the initial mass of the pure component oxides sample prior to the drop calorimetry process, then converted to a mole basis. Figure 2.11 shows the measured standard enthalpies of formation for both oxidized and reduced magnesium-manganese oxide. Figure 2.11. Acid Solution Calorimetry Results. Standard enthalpies of formation for the oxidized and reduced magnesium-manganese oxide phases measured from acid solution calorimetry for Mn/Mg molar ratio variants of 2/3, 1/1, and 2/1. The dissolved samples were previously quenched in the drop calorimeter to preserve the phase structure that exists at the drop temperature. Drop equilibrium conditions are specified in Figure 2.10 [31]. 1800 Oxidized Phase 1600 Reduced Phase 1400 1200 |ΔHf| (kJ mol-1) 1000 800 600 400 200 0 2/3 1/1 2/1 Mn/Mg From these values, the chemical energy stored is calculated per mole of oxygen released as ΔHf, RED -ΔHf, OX ΔHRxn = 2 y2 -y1 (2.11). 48 Here, y1 and y2 are the values of y at 1000 and 1500 °C, respectively. The results are then converted to a gravimetric basis based on the measured reactivities of each Mn/Mg ratio variation [29] with the conversion given by ΔHRxn &kJ mol-1 3 -1 -1 O2 '*VO2 &cmO2 gMgMnO '*1000 (g kg ) (2.12). 22400 (cm3O mol-1O2 ) 2 The values of y measured during the calorimetric process are shown in Figure 2.12. Figure 2.12. Excess Oxygen Content. Mass change parameters in the chemical formula for magnesium-manganese oxide measured during the calorimetric process. Recall that the definition of y is the excess oxygen present in the magnesium-manganese oxide material beyond the stoichiometric at a given equilibrium state [31]. The total energy density of the material is calculated by Eq. 2.1. The combined results for the sensible/latent, chemical reaction, and total energy densities of magnesium-manganese oxide for the thermochemical cycle operating between a reduction equilibrium of (TRED = 1500 °C, 49 PO2,RED = 0.2 atm) and an oxidation equilibrium of (TOX = 1000 °C, PO2,OX = 0.2 atm) are shown in Figure 2.13. Figure 2.13. Energy Density Results. Energy density measurements for the thermochemical cycle operating between equilibrium states of (TRED = 1500 °C, PO2,RED = 0.2 atm) and (TOX = 1000 °C, PO2,OX = 0.2 atm) for magnesium-manganese oxide with Mn/Mg molar ratios of 2/3, 1/1, and 2/1. Total energy density measurements are broken down into the sensible/latent and chemical reaction components [31]. The total energy density measurements obtained using this method are quite accurate, as the uncertainties in total energy density for this thermochemical cycle were found to be within 6.0%, 5.5%, and 6.1% of the measurement for 2/3, 1/1, and 2/1 Mn/Mg, respectively. 2.7. Optimizing for Maximum Volumetric Energy Density After completing the energy density measurements for the 2/3, 1/1, and 2/1 Mn/Mg variations, it is now important to find the optimal thermodynamic equilibrium states to bound the 50 thermochemical cycle that maximize the volumetric energy density of magnesium-manganese oxide. Maximizing on a volumetric basis is of interest because for a packed bed contained within a fixed volume, the amount of energy that can be stored in the system is limited by the volume of the reactor in which the bed is packed. Thus, energy densities in this section will be studied on a per unit volume basis. The CALPHAD predictions imply that some optimal mix of operating conditions and material composition exists. Figures 2.1 and 2.2 showed that increasing the magnesium content increases the available energy density below a maximum temperature of 1500 °C, likely because of a lowering of the onset of reduction temperature for higher Mg content materials, as the phase diagram of Panda et al indicates that the phase transition from spinel to monoxide occurs at lower temperatures with increasing magnesium [36]. Experimental measurements in Figure 2.13 confirmed that higher magnesium content leads to higher gravimetric energy density. However, the experimental reactivity data in Figure 2.1 show that increasing the magnesium content reduces the total amount of oxygen released, thereby reducing the available energy density. These two opposing effects imply the existence of an optimum Mn/Mg ratio that maximizes energy density. Further, the CALPHAD model predicts that a much larger reactivity and energy density can be achieved across all three Mn/Mg compositions by lowering the oxygen partial pressure during reduction. This suggests that the extent of reactivity at lower oxygen partial pressures of reduction should be investigated to find the maximum achievable energy density. In a practical reactor system, the extent to which the oxygen partial pressure can be reduced is governed by the pressure limits of large-scale industrial vacuum pumps. Common industrial scale vacuum pumps can achieve an absolute vacuum of 150 mbar (approx. 0.15 atm of total pressure) [53]. Thus, an operating oxygen partial pressure of 0.05 atm (approx. 0.25 atm of 51 total pressure for air) was chosen as an achievable starting point for reduction under low oxygen partial pressure. A variety of redox cycle conditions are referenced throughout this section. They are listed in Table 2.4 and will be referred to by number throughout the remainder of the section. Table 2.4. Redox Condition Reference Numbers. Reference numbers for redox conditions within Sec. 2.7. A “redox condition” refers to a given equilibrium (T,PO2) condition that is of interest for finding the optimal equilibrium bounds for a magnesium-manganese oxide thermochemical cycle. This chart is referenced throughout this section as “redox condition i" for simplicity, where “i” is an index referring to the condition number. Condition Number PO2, RED (atm) PO2, OX (atm) TRED (°C) TOX (°C) 1 0.2 0.2 1500 1000 2 0.2 0.2 1500 600 3 0.05 0.05 1500 600 4 0.05 0.2 1500 1000 5 0.01 0.01 1500 600 6 0.01 0.2 1500 1000 2.7.1. Experimental Methods 2.7.1.1. Variation of Molar Composition Reactive stability and reactivity for 2/3, 1/1, and 2/1 Mn/Mg under redox condition 1 have already been studied. To expand the range of studied compositions, reactive stability and reactivity for samples of Mn/Mg of 3/2 and 3/1 are studied and compared with 2/3, 1/1, and 2/1 Mn/Mg. Materials are prepared using the solid-state reaction method, where raw powders of the pure metal oxides are mixed according to the desired stoichiometric molar ratio, ball-milled for 24 hours to ensure even mixing, heat-treated at 1500 °C under air for 20 hours, then cooled to 1000 °C and heat treated for an additional 4 hours. The resulting solids are then crushed and 52 sieved to a particle size range of 125-180 μm. Samples of approx. 10 g (to ensure a representative sample) of the sieved particles are packed inside a ceramic alumina tube of approx. 25.4 mm OD and 19mm ID and cycled in a tube furnace for five stable cycles under redox condition 1. Cycling is performed at a rate of approx. 10 °C min-1, with 1 hour dwells at both TRED and TOX to ensure complete reduction/oxidation. Reaction rate data are acquired using mass-spectrometry to measure the outlet gas composition, converted to units of volume of oxygen reacted per gram of sample per minute, and then integrated by cycle, taking the average over five stable cycles as the total reactivity. 2.7.1.2. Lowering Oxygen Partial Pressure During Reduction The potential increase in reactivity by lowering the oxygen partial pressure during reduction is measured for 2/3, 1/1, 3/2, 2/1, and 3/1 Mn/Mg. The analytical approach for determining the increase in total energy density with increased reactivity follows from the illustration of a typical oxygen released versus temperature plot shown in Figure 2.14. 53 Figure 2.14. Illustration of Reactivity for Varied PO2. Illustration of the reactivity (VO2) of magnesium-manganese oxide during isobaric thermal reduction between some TOX and TRED =1500 °C, where TOX < TRED [32],[36], [37]. The analysis assumes that at some reference equilibrium temperature TOX,Ref, the phase composition of the material, and thus the excess oxygen content y, is the same for all PO2, RED. Measurements of VO2 from TOX,Ref to the maximum reduction temperature TRED under different PO2, RED facilitate the calculation of the increase in reactivity for a thermochemical cycle between (TRED, PO2, RED) and (TOX, PO2, OX) where PO2, RED < PO2, OX , relative to a cycle where PO2, RED = PO2, OX. Note that TRED > TOX > TOX,Ref . Samples of approximately 10 g are reacted according to redox condition 2, with 3-4 hour dwell times at both TRED and TOX to ensure complete reduction/oxidation. TOX,Ref was chosen as 600 °C to guarantee the same starting oxidation states [36], [37] and because of the nature of the tube furnace: 600 °C was achievable within a much faster cooling time than lower temperatures. Reaction rate data during the reduction step are collected using mass-spectrometry to measure the outlet gas composition and then integrated to calculate the reactivity for redox condition 2, defined as VO2, 2. The procedure is repeated using redox condition 3, from which VO2, 3 is 54 calculated. The increase in reactivity by lowering the oxygen partial pressure during reduction is given by, ∆VO2 = VO2, 3 – VO2, 2 (2.13). The reactivity for redox condition 4 is then calculated as, VO2, 4 = VO2, 1 + ∆VO2 (2.14) where VO2, 1 is the reactivity for redox condition 1. The same procedure is then repeated using redox conditions 1, 5, and 6 in order to determine the increase in reactivity for a cycle under redox condition 6, i.e., VO2, 6. Gravimetric energy densities per mass of fully oxidized sample are calculated from the measured reactivities and calorimetry data as 1000 g kg-1 ΔHS,i = ΔHRxn VO2,i 1 -1 2+ ΔHDrop, RED - ΔHDrop, OX (2.15) 22400 cm3 molO 2 where ΔHS,i is the energy density of the material operating under redox condition i, ΔHRxn is the energy of the redox reaction with units kJ molO2-1 calculated by Eq. 2.11 and E1. 2.12 (assuming ΔHRxn is independent of the extent of thermal reduction), VO2,i is the reactivity of the sample under redox condition i in cm3 g-1 calculated from the reactivity measurements, and ΔHDrop,RED and ΔHDrop,OX are the drop calorimetry measurements. Because the maximum volumetric energy density is of interest, ΔHS,i is converted to a volumetric basis using bulk density measurements of porous magnesium-manganese oxide samples. Samples of particle size 125-180 μm are formed into well-defined disks by heat treating at 1500 °C within alumina crucibles for 20 hours and subsequently oxidized at 1000 °C for 4 hours to form fully oxidized disks. Multiple diameter and thickness measurements were taken across each sample using a Vernier caliper, with the averages taken for calculations. Bulk 55 densities of the fully oxidized porous magnesium-manganese oxide disks are shown in Figure 2.15. Gravimetric energy densities are measured using the calorimetric procedure and combined with the density measurements to calculate volumetric energy densities for redox conditions 1, 4, and 6 for each Mn/Mg variation. Because of issues with the tube furnace at the time of this work, the sensible enthalpy differences for 3/2 and 3/1 Mn/Mg were estimated using the Neumann- Kopp procedure for mixed oxides described by Leitner et al. [54], neglecting dilation. Figure 2.15. Bulk Densities. Bulk densities of fully oxidized disks of magnesium-manganese oxide with Mn/Mg molar ratios of 2/3, 1/1, 3/2, 2/1, and 3/1. Disks were formed by sintering particles of 125-180 𝜇m in alumina crucibles at 1500 °C for 20 hours, followed by oxidation at 1000 °C for 4 hours. Densities are calculated from the mass and volume of each disk. These values are used for converting gravimetric energy density measurements to a volumetric basis [32] . 2.5 2 1.5 ! (g cm-3) 1 0.5 0 2/3 1/1 3/2 2/1 3/1 Mn/Mg 56 2.7.2. Results Reaction rate data acquired for five stable cycles under redox condition 1 for magnesium- manganese oxide with molar ratios of Mn/Mg of 3/2 and 3/1 are shown in Figure 2.16. Figure 2.16. Reactive Stability for New Mn/Mg Samples. Reaction rates over five stable cycles of magnesium-manganese oxide of particle size 125-180 μm under redox condition 1 for molar ratios of Mn/Mg of (a) 3/2 and (b) 3/1 [32]. Samples appear to exhibit cyclically stable reactive behavior over five cycles. The total reactivity is found by integrating each cycle over time, then taking an average of the reactivity per cycle. Reactivities were found to be 35.0 ± 3.4 and 18.3 ± 3.8 cm3 g-1 for 3/2 and 3/1, respectively [32]. (a) 57 Figure 2.16 (cont’d). (b) The total reactivity averaged over the five cycles for each sample was found to be 18.3 ± 3.8 and 35.0 ± 3.4 cm3 g-1 for 3/1 and 3/2, respectively. The increases in reactivity and subsequent increases in volumetric energy density achievable by lowering the oxygen partial pressure during reduction for Mn/Mg samples of 2/3, 1/1, 3/2, 2/1, and 3/1 are shown in Figures 2.17 and 2.18, respectively. 58 Figure 2.17. Measured Reactivity for Varied PO2. Reactivities of magnesium-manganese oxide with Mn/Mg molar ratios of 2/3, 1/1, 3/2, 2/1, and 3/1 for a thermochemical cycle with an oxidation equilibrium state of (TOX = 1000 °C, PO2,OX = 0.2 atm) and reduction equilibrium states of (TRED = 1500 °C, PO2,RED) where PO2,RED = 0.2, 0.05, or 0.01 atm. Reactivity increases significantly for all samples with decreasing PO2,RED [29], [31], [32]. 70 0.2 atm 0.05 atm 0.01 atm 60 50 40 VO2 (cm3 g-1) 30 20 10 0 2/3 1/1 3/2 2/1 3/1 Mn/Mg 59 Figure 2.18. Volumetric Energy Densities at Varied PO2. Volumetric energy densities of magnesium-manganese oxide with Mn/Mg molar ratios of 2/3, 1/1, 3/2, 2/1, and 3/1 for a thermochemical cycle with an oxidation equilibrium state of (TOX = 1000 °C, PO2,OX = 0.2 atm) and reduction equilibrium states of (TRED = 1500 °C, PO2,RED) where PO2,RED = 0.2, 0.05, or 0.01 atm. Volumetric energy density increases significantly for all samples with decreasing PO2,RED [31], [32] . 3000 0.2 atm 0.05 atm 0.01 atm 2500 2000 ΔHS (MJ m-3) 1500 1000 500 0 2/3 1/1 3/2 2/1 3/1 Mn/Mg Figure 2.17 and 2.18 demonstrate that increasing the manganese content of the material above a 2/1 Mn/Mg ratio yields a greatly diminished reactivity under redox condition 1, resulting in a lower energy storage density for 3/1 Mn/Mg. Particularly, 3/1 Mn/Mg exhibits a 20.5% reduction in volumetric energy density relative to 2/1 Mn/Mg. This indicates that materials with a high manganese content (> 2/1 Mn/Mg) are not practical for engineering applications due to substantially lower energy density. It is also clear from the bulk density measurements and the 60 figures that the magnesium content and the physical density have the biggest effect on volumetric energy density. Increasing magnesium content increases the energy stored as sensible heat, while conversely decreasing the physical density. For 2/3, 1/1, and 2/1 Mn/Mg, the differences in reactivity and standard enthalpy of formation are small enough that the volumetric energy densities are strongly dependent on the bulk material density. Thus, for Mn/Mg ratios between 2/3 & 2/1 Mn/Mg, the optimal volumetric energy density seems to depend on the balance between magnesium content and bulk density, appearing to peak at 1/1 Mn/Mg and indicating that the optimal Mn/Mg ratio for maximum volumetric energy density lies somewhere around this composition. Figures 2.17 and 2.18 also demonstrate the gains in reactivity and energy density achievable by lowering the oxygen partial pressure to both 0.05 and 0.01 atm during reduction. Results indicate that this is an effective approach for increasing the energy density of magnesium-manganese oxide, as energy densities increased noticeably for each material composition. The effect is strongest for high manganese content materials: reducing under PO2,RED = 0.05 atm improves the energy density of 3/1 Mn/Mg by approximately 42% from that measured for PO2,RED = 0.2 atm, indicating that Mn/Mg ratios above 2/1 are substantially more viable for low oxygen partial pressure reduction conditions. The trend continues for all Mn/Mg samples when lowering PO2,RED further to 0.01 atm. Again, we observe significant increases in energy densities, with the greatest gains in energy density observed for high Mn-content materials. A point of particular interest is that volumetric energy densities for 1/1, 3/2, and 2/1 Mn/Mg all approach approximately 2300 MJ m-3, indicating that when operating under very low PO2,RED, volumetric energy densities will reach a maximum and thus the optimal Mn/Mg ratio could be selected from amongst 1/1, 3/2, and 2/1 Mn/Mg based on other design factors, such as 61 the fastest reaction kinetics. However, 1/1 Mn/Mg is consistently the optimal choice for maximizing volumetric energy density under all PO2,RED and is therefore the optimal Mn/Mg ratio for a magnesium-manganese oxide TCES reactor. Figure 2.15 and 2.16 demonstrate that increasing the manganese content of the material above a 2/1 Mn/Mg ratio yields a greatly diminished reactivity under redox condition 1, resulting in a lower energy storage density for 3/1 Mn/Mg. Particularly, 3/1 Mn/Mg exhibits a 20.5% reduction in volumetric energy density relative to 2/1 Mn/Mg. This indicates that materials with a high manganese content (> 2/1 Mn/Mg) are not practical for engineering applications due to substantially lower energy density. It is also clear from the bulk density measurements and the figures that the magnesium content and the physical density have the biggest effect on volumetric energy density. Increasing magnesium content increases the energy stored as sensible heat, while conversely decreasing the physical density. For 2/3, 1/1, and 2/1 Mn/Mg, the differences in reactivity and standard enthalpy of formation are small enough that the volumetric energy densities are strongly dependent on the bulk material density. Thus, for Mn/Mg ratios between 2/3 & 2/1 Mn/Mg, the optimal volumetric energy density seems to depend on the balance between magnesium content and bulk density, appearing to peak at 1/1 Mn/Mg and indicating that the optimal Mn/Mg ratio for maximum volumetric energy density lies somewhere around this composition. Figures 2.15 and 2.16 also demonstrate the gains in reactivity and energy density achievable by lowering the oxygen partial pressure to 0.05 and 0.01 atm during reduction. Results indicate that this is an effective approach for increasing the energy density of magnesium-manganese oxide, as energy densities increased noticeably for each material composition. The effect is strongest for high manganese content materials: reducing under 62 PO2,RED = 0.05 atm improves the energy density of 3/1 Mn/Mg by approximately 42% from that measured for PO2,RED = 0.2 atm, indicating that Mn/Mg ratios above 2/1 are substantially more viable for low oxygen partial pressure reduction conditions. The observed trends continue for all Mn/Mg samples if lowering PO2,RED further to 0.01 atm. Again, we observe significant increases in energy densities, with the greatest gains in energy density observed for high Mn-content materials. A point of particular interest is that volumetric energy densities for 1/1, 3/2, and 2/1 Mn/Mg all approach approximately 2300 MJ m- 3 , indicating that when operating under very low PO2,RED, volumetric energy densities appear to reach a maximum and thus the optimal Mn/Mg ratio could be selected from amongst 1/1, 3/2, and 2/1 Mn/Mg based on other design factors, such as the fastest reaction kinetics. However, 1/1 Mn/Mg consistently appears to be the optimal choice for maximizing volumetric energy density under all PO2,RED and is therefore the optimal Mn/Mg ratio for a magnesium-manganese oxide TCES reactor. 2.8. Comparison of Energy Densities with CALPHAD Predictions Recall that in Sec. 2.2, the Panda et al [36] CALPHAD model for the Mg-Mn-O system is used to predict the reactivities and energy densities of 2/3, 1/1, and 2/1 Mn/Mg for reduction at 1500 °C at a few oxygen partial pressures of reduction, with the reference oxidation equilibrium state of (TOX = 1000°C, PO2,OX = 0.2 atm). Now that these data have been measured, it is of interest to compare the data with the CALPHAD predictions. Figures 2.19 and 2.20 give these comparisons for the reactivity and energy density results, respectively. Note that the measured energy densities for a reduction equilibrium state of (TRED = 1500°C, PO2,RED = 0.01 atm) are converted from the volumetric basis to a gravimetric basis using the reported bulk densities of the 125-180 𝜇m samples. 63 Figure 2.19. Comparison with CALPHAD, Reactivities. Comparison of the measured reactivities for magnesium-manganese oxide with Mn/Mg molar ratios of 2/3, 1/1, and 2/1 with CALPHAD predictions for a thermochemical cycle with an oxidation equilibrium state of (TOX = 1000 °C, PO2,OX = 0.2 atm) and reduction equilibrium states of (TRED = 1500 °C, PO2,RED) where PO2,RED = 0.2 or 0.01 atm [29], [31], [32], [36]. When compared to the measured values, the CALPHAD predictions become increasingly less accurate as the composition of the material moves further away from the stoichiometric MgMn2O4 on which the model predictions are based. 70 CALPHAD, 0.2 atm CALPHAD, 0.01 atm 60 Measured, 0.2 atm Measured, 0.01 atm 50 40 VO2 (cm3 g-1) 30 20 10 0 2/3 1/1 2/1 Mn/Mg 64 Figure 2.20. Comparison with CALPHAD, Energy Densities. Comparison of the measured gravimetric energy densities for magnesium-manganese oxide with Mn/Mg molar ratios of 2/3, 1/1, and 2/1 with CALPHAD predictions for a thermochemical cycle with an oxidation equilibrium state of (TOX = 1000 °C, PO2,OX = 0.2 atm) and reduction equilibrium states of (TRED = 1500 °C, PO2,RED) where PO2,RED = 0.2 or 0.01 atm [29], [31], [32], [36]. When compared to the measured values, the CALPHAD predictions become increasingly less accurate as the composition of the material moves further away from the stoichiometric MgMn2O4 on which the model predictions are based. 1800 CALPHAD, 0.2 atm CALPHAD, 0.01 atm 1600 Measured, 0.2 atm Measured, 0.01 atm 1400 1200 ΔHS (kJ kg-1) 1000 800 600 400 200 0 2/3 1/1 2/1 Mn/Mg The figures show that the CALPHAD model seems to perform reasonably well for 2/1 Mn/Mg, as the predicted reactivities and energy densities are consistently within the experimental margin of error. However, this begins to break down with increasing magnesium content. For 1/1 Mn/Mg, the model was within the margin of error for the energy density 65 measurement for a reduction equilibrium state of (TRED = 1500°C, PO2,RED = 0.2 atm), but significantly outside the margin of error for the remaining measurements. For 2/3 Mn/Mg, the model was quite close to the experimental reactivity for a reduction equilibrium state of (TRED = 1500°C, PO2,RED = 0.2 atm), but is consistently low for all remaining 2/3 Mn/Mg measurements. The interesting trend is that the CALPHAD model seems capable of reasonably reproducing the Mg-Mn-O system in these high temperature ranges when the composition most closely resembles the stoichiometric spinel of MgMn2O4, but generally underestimates the amount of oxygen the material can absorb with higher Mg-content, thus underpredicting the energy density that the system is capable of when more Mg is added. This might have something to do with the non-stoichiometric effects [29] and warrants a further comprehensive thermodynamic study of the Mg-Mn-O system at temperatures above 1000 °C to determine the cause and develop a more accurate thermodynamic model. 2.9. Bulk Oxidation Reaction Kinetics Up to this point, only the total energy density and the optimal thermodynamic equilibrium conditions for maximizing volumetric energy density have been identified. However, this lacks any quantified information as to the rate at which the material releases that energy during the oxidation reaction. To model the energy discharge rate of a reacting packed bed of magnesium-manganese oxide particles, the oxidation rate needs to be quantified. In most TCES redox metal oxide studies, kinetics are analyzed using very small samples (order of milligrams) in thermogravimetric analyzers [8]-[28]. However, for large-scale engineering design purposes it may be beneficial to have information about the oxidation kinetics of bulk samples. This section summarizes the work done to quantify the oxidation kinetics of bulk samples of magnesium- manganese oxide particles in the temperature range 1000-1500 °C. 66 2.9.1. Experimental Setup and Procedure To examine a few different molar compositions, samples of 2/3, 1/1, and 2/1 Mn/Mg are synthesized using the solid-state reaction approach. Raw component powders are pre-mixed in the given molar ratios and then milled for 24 hours in a porcelain ball-mill jar with zirconia grinding media. The mixed raw powders are then heat-treated in alumina crucibles at 1500 °C for an additional 24 hours. The resulting sintered solids are then crushed, ground, and sieved to the particle size range of 125-180 𝜇m. Approx. 10 g of sieved particles are packed in the center of an alumina ceramic tube (OD 25.4 mm, ID 19mm) and supported with Buster M-35 high- temperature insulation from Zircar Zirconia, which also reduces end losses from the sample. Sample temperature at the inlet is recorded using an alumina-sheathed B-type thermocouple inserted within the tube. The composition of the outlet gas mixture is recorded with a Hiden HPR-20 quadrople mass-spectrometer using a Faraday scanning technique. Samples are heated to 1500 °C in a Sentrotech STT-1700°C-6 tube furnace under a constant 1 SLPM flow of 99.999% pure argon and held for 4 hours to ensure full thermal reduction, then cooled to the desired oxidation temperature while maintaining the argon flow. The sample is then held for 1 additional hour to ensure thermal equilibrium. At this point the gas flow mixture is changed to consist of 80 SCCM of argon and 20 SCCM of oxygen such that the sample is exposed to a sudden oxygen partial pressure of 0.2 atm, approximately the equivalent of air. The mass-spec monitors the oxygen concentration in the outlet gas until the outlet composition matches that of the inlet, indicating that the reaction has completed. A previous series of blank runs using a similar mass of alumina particles packed under the same conditions showed that without any reaction, the oxygen signal at the mass-spec jumped from 0% to 20% oxygen concentration almost immediately after adjusting the flow rates, implying that any delay 67 in oxygen detection by the mass-spec during the experiments is entirely due to oxidation of the sample. The blank runs thus serve as a baseline for the oxidation runs, where the difference between the two gives the oxidation rate and the integrated rate gives the total oxygen consumed during the run. This process is done isothermally at increments of 100 °C between 1000-1500 °C to adequately cover the thermochemical cycling range. A schematic of the setup is shown in Figure 2.21. Figure 2.21. Oxidation Kinetic Experiments, Setup. Experimental setup for the bulk oxidation kinetic experiments [33]. 2.9.2. Kinetic Modeling and Data Extraction The modeling work for extracting the kinetic rate law and parameters is briefly summarized here to include the results while appropriately reflecting the author’s involvement in this portion of the work. It is important to note that all modeling and data extraction work for the 68 bulk oxidation kinetics was performed by Dr. Kelvin Randhir and is outlined in much more detail in [33]. The kinetic rate law is derived using Schmalzried’s theory of internal oxidation [55]-[58] with a highly-simplified stochiometric version of the magnesium-manganese oxide reaction system and assuming the spherical particle analogy in which each particle oxidizes under the same conditions. Defining the extent of reaction 𝛼 as the ratio of the mass of oxygen absorbed at time t to the total final mass of oxygen absorbed and applying the above theory results in the kinetic rate law, given by Eq. 2.15. The full derivation can be found in [33]. Ea 2 PO - 2, Ko e RT (1- α)3 + dα PRef dt = 1 (2.15). (1- (1-α)3 ) + S Here, Ko, Ea, and S are the three kinetic parameters to be fit. To extract them from the experimental data, a 1D plug-flow conservation of species model is derived following the approach of Mehdizadeh et al [59] to model oxygen concentration for the bulk system, as oxygen is the only reacting species. The system of partial differential equations is given by Eq. 2.16 and 2.17, with the full derivation available in [33]. ∂γO ∂γO RTṙ O 2 2 ∂t = -u ∂z + ϵPM 2 (1-γO ) (2.16), O2 2 dα ṙ O2 = -mO2 ,∞ dt (2.17). Here, T is the temperature of the gas, R is the universal gas constant, P is the pressure, MO2 is the molar mass of oxygen, 𝜖 is the bulk porosity of the sample, 𝛼 is the extent of reaction, 𝛾O2 is the mole fraction of oxygen, u is the interstitial axially velocity of gas, ṙ O2 is the rate of oxygen consumption, and mO2 ,∞ is the maximum extent of oxygen absorbed by the material at a fixed temperature T. The model is discretized using a 1D finite difference upwind scheme and solved 69 explicitly over time for the oxygen concentration at the outlet node, starting at a fully reduced state and proceeding to a fully oxidized state. A global least-squares-fit optimization of the three kinetic parameters is done using fmincon in MATLAB to optimize the rate parameters simultaneously across all datasets. 2.9.3. Results The optimized kinetic parameters are given in Table 2.5. The total oxygen absorbed for each material sample during each isothermal run is shown in Figure 2.22. The experimental and modeled outlet oxygen flow data for each sample are shown in Figure 2.23. Table 2.5. Oxidation Kinetic Parameters. Optimized bulk oxidation kinetic parameters for magnesium-manganese oxide between 1000-1500 °C [33]. Ko is the rate constant, Ea is the activation energy, and S is a constant parameter in the simplified rate law given by Eq. 2.15. Kinetic Parameter Units 2/3 Mn/Mg 1/1 Mn/Mg 2/1 Mn/Mg Ko s-1 0.866 0.852 0.641 Ea kJ mol-1 57.5 57.1 56.2 S - 0.0474 0.0587 0.0512 70 Figure 2.22. Oxidation Kinetic Experiments, Total Absorbed. Total oxygen absorbed by each magnesium-manganese oxide sample during the isothermal oxidation runs [33]. Values are found by integrating the reactivity rate data acquired by the mass-spec over the course of each run. The y-axis has units of 𝜇mol of oxygen gas per gram of solid material. 3500 2/3 Mn/Mg 3000 1/1 Mn/Mg 2500 2/1 Mn/Mg mO2,∞ (!molO2 g-1) 2000 1500 1000 500 0 900 1000 1100 1200 1300 1400 1500 1600 T (°C) 71 Figure 2.23. Oxidation Kinetic Experiments, Data and Model Fits. Experimental and modeled oxygen outlet flow for all isothermal runs for magnesium-manganese oxide with Mn/Mg molar ratio of (a) 2/3, (b) 1/1, and (c) 2/1 [33]. The model results using the globally- optimized kinetic parameters fit reasonably well across the 1000-1500 °C temperature range for all samples, with the exception of that at 1400 °C. It is hypothesized that this deviation is due to a phase change during the oxidation run, as all three samples are near the phase boundary line in the phase diagram of Panda et al [36]. However, validation of this requires much further study that is beyond the scope of this work. (a) 72 Figure 2.23 (cont’d). (b) 73 Figure 2.23 (cont’d). (c) It is clear from Figure 2.22 that the bulk kinetic model is incapable of accurately modeling oxidation at 1400 °C. This is possibly because of an active phase transition around 1400 °C: the Mg-Mn-O phase diagram of Panda et al [36] shows that these Mn/Mg compositions are near the spinel-to-monoxide phase transition around this temperature. However, no conclusions should be drawn, and this apparent deviation in the oxidation mechanism should be 74 thoroughly investigated. Overall, the average R2 value across all curve-fits for each sample was 0.985, 0.976, and 0.970 for 2/3, 1/1, and 2/1 Mn/Mg respectively. It is concluded that while further investigation of the oxidation mechanisms at 1400 °C is necessary, the bulk oxidation model is suitably accurate (R2 > 0.95) to model the oxidation component of the energy discharge process in a reacting bed of magnesium-manganese oxide between 1000-1500 °C. 2.10. Summary This chapter determined the relevant thermochemical properties of magnesium- manganese oxides for a TCES reactor operating between 1000-1500 °C. Energy densities for 2/3, 1/1, and 2/1 Mn/Mg materials were predicted using CALPHAD modeling. A calorimetric method combining drop and acid solution calorimetry was developed, validated, and used to measure the energy density of the same three materials assuming oxidation at 1000 °C, reduction at 1500 °C, and a constant oxygen partial pressure of 0.2 atm. In light of CALPHAD predictions, an investigation was done to determine the maximum extent of energy density achievable by reducing the partial pressure of oxygen during the reduction step. Compositions of 3/2 and 3/1 Mn/Mg were also studied to expand the compositional range. The maximum achievable volumetric energy density for magnesium-manganese oxide is possible by selecting the molar composition of Mn/Mg 1/1 for a thermochemical cycle operating between the reduction equilibrium state of (TRED = 1500 °C, PO2,RED = 0.01 atm) and the oxidation equilibrium state of (TOX = 1000 °C, PO2,OX = 0.2 atm). The 1/1 Mn/Mg variant is selected because it has the highest volumetric energy density, making it the optimal choice for a fixed volume reactor vessel. Finally, a bulk oxidation kinetic model was derived based on the theory of internal oxidation using a 1D plug flow species conservation approach. The resulting global model had average R2 75 values > 0.95 and is suitable for coupling with the packed bed heat transfer effects for modeling a non-isothermal reacting packed bed of magnesium-manganese oxide. 76 CHAPTER 3: OVERVIEW OF POROUS MEDIUM HEAT AND MASS TRANSFER THEORY The thermodynamic measurements outlined in Chapter 2 provide most of the information needed to model a reactive packed bed of magnesium-manganese oxide particles. The total energy density, the optimal thermodynamic equilibrium states to bound the thermochemical cycle within which the reactor should operate, and the rates at which the material releases the chemical reaction energy were all quantified. However, a model for the transport of that energy from the material within the bed to the bed outlet and on to the power generation system outside the bed is still needed. To model this portion of the energy discharge step, the heat and mass transport equations that govern the system need developed. As previously stated, an ideal magnesium-manganese oxide TCES reactor will consist of a packed bed of particles to maximize the surface area that can interact with the flow of air (given that air is both the reactant and the working fluid). Thus, the transport model must come from the heat and mass transport theory for flow through a porous medium. The remainder of this chapter overviews the aspects of porous medium heat and mass transport theory that are relevant to flow through a packed bed. A porous medium can be thought of as a composite material with two phases. The first is a matrix-type phase consisting of the solid material. The second is an inclusion-type phase where the “inclusions” are the empty spaces within the solid matrix, called “voids” or “pores”. The pores are filled by the fluid in which the solid is immersed, such as water or air. When the pores are interconnected, the fluid is free to flow along paths through the pore space around the solid matrix. Porous media are often used by engineers when designing systems involving heat transfer and/or chemical reaction because of the significant gains in the available surface area relative to a non-porous material. However, this enhanced surface area leads to complex geometries that pose an obstacle to the theoretical treatment of flow through the pores. Each 77 porous medium often has a unique geometric configuration, making a unified pore-level theory a difficult and as yet unaccomplished task. To solve this problem, scientists and engineers have developed a volume-averaging approach that allows the solid and fluid phases to be modeled as two separate continua in the same domain. This approach has produced a collection of theoretical equations that allows for the average behavior of fluid flow and heat transfer within the porous medium to be modeled at the length scales relevant to the macroscopic domain. This chapter outlines the portions of porous media heat and mass transport theory relevant to modeling the energy discharge from a reactive bed of magnesium-manganese oxide particles. The volume-averaging theory is summarized and condensed into its core elements. The resulting fluid mechanics, conductive heat transfer, and convective heat transfer averaged on the macroscopic-scale for single-phase flow through a homogenous porous medium are summarized, and an approach for modeling the heat and mass transfer within the packed bed is chosen. 3.1. Volume Averaging Theory The volume averaging theory of porous media transport is derived in Whitaker’s The Method of Volume Averaging [60]. In the text, he describes the theory by way of application to the species, momentum, and energy conservation equations. He derives the macro-scale volume- averaged conservation equations step by step, starting from the base forms at the pore level of a two-phase porous medium. Because the derivation is quite in-depth, only the important elements are summarized in this work. To begin, a length scale is chosen at which a reasonable average of a given property can be achieved to form the basis for a representative elementary volume (REV). An REV is defined as the smallest volume over which a measurement can be made and still give a representative measurement of some quantity within the domain. The REV is selected such that the local 78 fluctuations of the state variables and system properties at a given point average out over the REV to a mean value that is measurable at the length scales of interest. Every porous medium has a different length scale that forms the REV and must be chosen with care. An illustration of a hypothetical REV within a porous domain is shown in Figure 3.1. Figure 3.1. Representative Elementary Volume. Hypothetical representative elementary volume over which the continuum equations are averaged within a porous domain. The length scale of the REV is selected such that l << L, yet the portion of the porous medium contained within the REV gives a reasonable average representation of the local geometry. Application of the volume averaging laws to the Navier-Stokes equations on the scale of the REV will yield an average value of a given property within the REV shown. After defining a suitable REV, the pore-level conservation equations are spatially smoothed over the REV. This is done using the volume averaging laws outlined below. First, the phase fraction is defined as V 𝜖. = Vi (3.1) where Vi is the volume of the constituent i within the REV, V is the total volume of the REV, and 𝜖. is defined as the fraction of the total REV that consists of constituent i. The phase fractions of all phases must sum to 1. For single phase flow through a homogeneous porous medium, there are only two separate phase fractions and thus they can be defined as 𝜖 and 1-𝜖. 79 Next, one must define the two types of averages: superficial and intrinsic. The superficial average is defined as the average of a property ϕi over the entire REV domain volume V. In mathematical terms, this is given as 〈ϕi 〉 = ∫V ϕi dV (3.2). Here, the angled brackets denote that the enclosed quantity is a volume-averaged value. The intrinsic average is defined as the average of a property ϕi over only the respective phase volume Vi. Mathematically this is given as 〈ϕi 〉i = ∫V ϕi dVi (3.3). i The intrinsic and superficial averages are related through the relationship, 〈ϕi 〉 = ϵi 〈ϕi 〉i (3.4). Eq. 3.4 means that the superficial average of phase i is equal to the intrinsic average of phase i redistributed over the REV domain by weighing by the volume fraction of phase i. The superficial and intrinsic averages are more easily understood with a visualization. Figure 3.2 illustrates the concept of intrinsic and superficial averages for the flow of a fluid throw the REV shown in Figure 3.1. 80 Figure 3.2. Intrinsic and Superficial Averages. Illustration of the intrinsic and superficial averages from Eq. 3.2 and 3.3 as applied to a fluid flow through the representative elementary volume shown in Figure 3.1. The leftmost REV shows the fluid flow at the pore level governed by the Navier-Stokes equations. The middle REV represents the intrinsic average, where the pore-level flow velocity vectors are averaged over the highlighted volume of the fluid phase. The rightmost REV represents the superficial average, where the pore-level flow velocity vectors are averaged over the entire REV. In addition to the definitions of the phase fraction, the superficial average, and the intrinsic average, one must also make use of the spatial averaging theorem, given by 1 〈∇ϕi 〉 = ∇〈ϕi 〉 + ∫ nij ϕi dA (3.5). V Here, nij is the normal vector to the boundary surface between phase i and phase j, and A is the total surface area along the boundary between phase i and phase j. In words, the spatial averaging theorem says that the average of the gradient of some property ϕi can be represented as the sum of the gradient of the average property over the REV and the surface integral of the property ϕi along the phase boundary divided by the REV’s total volume. This allows terms to be reconfigured over the course of the volume-averaging treatment into forms that are more easily understood, observed, and measured. The treatment also makes use of the spatial decomposition rule, given by ϕi = 〈ϕi 〉i + 8 ϕi (3.6). 81 This expression allows the property ϕi to be broken down into the intrinsic average 〈ϕi 〉i and the local deviation from the intrinsic average, 8 ϕi . This is important because it separates the length scales of the average property (on the order of the REV) and the deviations (order much less than the REV), allowing some terms involving the deviations from the average to be neglected on the basis of negligible length scale. Whitaker [60] notes that every volume averaging problem has to handle three terms that appear from application of the spatial averaging laws: (1) volume averaged properties inside phase boundary surface integrals, (2) properties averaged over the phase boundary surface instead of over the REV, and (3) the presence of local spatial deviations in the equations. Each problem is handled in a different manner using a combination of Taylor series expansions, length scale order of magnitude comparisons, and some geometric theorems developed by Quintard and Whitaker [61]. Some boundary surface integral terms come to represent the cross-phase exchange terms (interstitial convective heat transfer, species transport between phases, etc.). Finally, the last spatial deviations are removed by solving the closure problem, making use of all the volume- averaging laws and length scale order of magnitude comparisons used in the earlier steps. The result is a collection of continuum equations at the macroscopic scale of the REV capable of modeling the heat and mass transport phenomena occurring within the porous medium at observable scales. Whitaker applies the treatment to the energy equation assuming local thermal equilibrium between the two phases [61], eliminating interstitial convection from the equation. In a separate work, Quintard and Whitaker perform the same treatment for the case of two distinct phases in non-local thermal equilibrium, resulting in a generalized form of the two-phase energy model containing interstitial convection [62],[63]. This work is interested in the non-local thermal 82 equilibrium equations because they allow for interstitial convective heat transfer between the two phases. It is important to have a general knowledge of where the volume-averaged equations come from if one is to use them. For more complex systems with heterogeneous geometry or multi-phase flow, it may be necessary to derive a unique set of equations using the volume- averaging process. However, for the purposes of this study, single-phase flow through homogeneous porous media under non-local thermal equilibrium is well-studied with established simplified volume-averaged governing equations, thus making a unique derivation unnecessary. Here, homogeneous means that the porous medium is isotropic (no directional-dependent properties) and disordered (no channeling or changes of the average flow behavior due to ordered structuring of the matrix). The remaining sections of this chapter summarize the fluid mechanics and energy transport relative to this particular geometric case. 3.2. Fluid Mechanics in Homogenous Porous Media with Single-Phase Flow The fundamental law for fluid mechanics in porous media is the Darcy law [64],[65] , developed by Henry Darcy in 1856 from his work with the flow of water through packed columns of sand (a disordered, isotropic porous medium). This law is given mathematically as dP μ - dz = K uD (3.7) dP where dz is the pressure drop along the column, 𝜇 is the dynamic viscosity of the fluid flowing through the column, K is the permeability of the porous medium (a measure of the ease with which a fluid can flow through the geometry of a given porous medium with units m2 s-1), and uD is the Darcy velocity (i.e., the superficial velocity) of the fluid. In words, Darcy’s law says that the pressure drop through a porous medium is proportional to the superficial velocity. For isotropic, disordered media, permeability is a constant. Non-isotropic media must treat 83 permeability as a tensor with different magnitudes depending on the direction of flow. Different theoretical models exist to predict permeability for a given porous medium [64], but are beyond the scope of this discussion. As the flow increases, the behavior deviates from Darcy’s law. This is attributable to increased inertial or non-linear effects within the pores. For flow through porous media, the appropriate length scale for the Reynolds number is the average pore diameter (equivalent to the average particle diameter in a packed bed). When defined as such, the flow behavior begins to deviate from Darcy’s law at approx. a Reynolds number of 1 [66]. Beyond this point, inertial/non- linear effects become increasingly important until becoming the dominant force at high Darcy velocities. Kaviany categorizes the four flow regimes at the pore level as follows: [64] 1) Redp < 1, Darcy regime, Viscous forces dominate, linear laminar flow 2) 1-10 < Redp < 150, Inertial/Non-linear regime, inertia/non-linear effects increasingly dominates with increasing flow, non-linear laminar flow 3) 150 < Redp < 300, Transition to pore-level turbulence 4) 300 < Redp , Fully turbulent at the pore-level. Clearly, Darcy’s law loses modeling relevance as the flow rate increases. Correction terms developed by Brinkmann and Forcheimer [64],[66], [67] were added to Darcy’s law to account for the inertial effects. Ergun [68] took a similar approach, using hydraulic radius theory to develop a relationship that accounts for these effects and can accurately quantify the pressure drop within packed beds at a range of flow rates. The most versatile of the Darcy law extensions is the Darcy-Brinkmann-Forcheimer generalized flow model, where the particular form shown in Eq. 3.8 was derived from the volume-averaging treatment by Hsu and Cheng [66], [69]. For an isotropic, disordered porous medium with all properties appropriately volume-averaged, the equation is 84 〈ρ〉 ∂〈u1 〉 (〈u1〉*∇)〈u1〉 μ 〈ρ〉 b 〈u1〉 |〈u1 〉| 9 + : = 〈G 〉 - ∇〈P〉 + μ∇2 〈u"〉 - K 〈u" 〉+ (3.8). ϵ ∂t ϵ K0.5 Here, 𝜌 is the fluid density, 𝜖 is the porosity of the porous medium (i.e., the volume fraction of the fluid phase), G represents body forces, 𝑢" is the velocity vector, and b the Forcheimer friction factor. The Darcy-Brinkmann-Forcheimer equation is a general macroscopic flow model that can be used when one needs to model flow through an isotropic, disordered porous medium beyond the Darcy regime. However, if the medium of interest is sufficiently dense such that solving Eq. 3.8 yields a plug-flow profile, then the system can be modeled as 1D where the Darcy velocity uD is sufficient and can be calculated from the mass flux through the system. It is also notable that while turbulence certainly occurs in the flow at the pore-level, dense porous media do not develop turbulence at the macroscopic level at which Eq. 3.8 is derived because the dense solid matrix structure prevents macroscopic turbulent eddies from forming [70],[71] . Thus, this model is more than sufficient for a wide range of flow rates and turbulence modeling only needs to be considered when modeling the fluid flow directly at the pore scale [70],[71] . 3.3. Conduction Heat Transfer in Homogeneous Porous Media Heat conduction at the REV scale within a porous medium is governed by an effective thermal conductivity, keff, in which the rates of heat conduction through the two phases are averaged over the REV. Effective thermal conductivities are unique to a given material and depend on four factors [64]: (1) the ratio of the solid-to-fluid conductivities ks/kf, (2) the extent of continuity of the solid matrix (a consolidated structure, a tightly packed bed, a loosely packed bed, etc.), (3) the contact resistance between particles if a non-consolidated medium, and (4) the magnitude of the Knudsen number for the gas phase (if large, the bulk gas conductivity cannot be used in a fluid thermal equation). 85 Effective thermal conductivities are typically measured for the porous medium of interest. However, theoretical models do exist. The most general theory applies the volume- averaging treatment to the local thermal equilibrium single-equation energy model to derive parameters for the effective thermal conductivity tensor. Further, a collection of theoretical and experimental studies have been performed to derive general effective thermal conductivity models [64], review of which is beyond the scope of this work. Kaviany gives a summary of the collection of models for beds of spheres and suggests Hadley’s Weighted Average model as a suitable choice [64], [72]. This model is given by the following collection of equations [64], [72] : f0 = 0.8+0.1ϵ (3.9) log(α0 ) = -1.084 - 6.778(ϵ - 0.298), 0.298 ≤ ϵ ≤0.580 (3.10) 2 ks k (1 + 2ϵ)ks keff ϵf0 + kf !1 - ϵf0 # 2+ks , (1 - ϵ)+ kf kf = (1 - α0 ) ks + α0 f (2 + ϵ) ks (3.11). 1- ϵ 41 - f0 5 + +1-ϵ kf ϵ !1 - f0 # kf Here, keff is the effective thermal conductivity of the bed, kf is the thermal conductivity of the fluid phase, ks is the thermal conductivity of the solid phase, 𝜖 is the bulk porosity of the porous medium, and f0 and 𝛼0 are model-specific parameters. This model will be used in Chapter 4 when working with packed beds of spheres. Effective thermal conductivities for beds of magnesium- manganese oxide particles have been measured in a separate study [30] and will be used in Chapter 4 when working with packed beds of magnesium-manganese oxide. 3.4 Convective Heat Transfer in Homogeneous Porous Media Convective heat transfer within a porous medium is yet another complex process. The rate of convection from the solid matrix to the fluid flow through the pores depends on flow rate, geometry of the medium (porosity, pore size), and the thermophysical properties of both phases. The volume-averaging treatment of Whitaker [60] focuses on a single-medium energy equation in 86 which local thermal equilibrium is assumed, eliminating the convective heat transfer from the equation. The much longer treatment by Quintard and Whitaker [62] removes the local thermal equilibrium assumption, deriving the general two-equation energy model for a two-phase porous medium. A follow-up work numerically determines the effective properties for certain specified geometries [63]. These results apply only to the specific geometries for which they were computed; therefore, most studies utilize highly simplified forms of the volume averaged equations that are sufficient for modeling the energy transport within homogeneous porous media. There are three such models of increasing complexity in which the two phases are treated as separate and overlapping continua coupled through interstitial convection. The first is the Schumann model [64],[73] in which the temperature distribution in the porous medium is modeled assuming a uniform plug flow through the bed, negligible radial temperature gradients, uniform temperature for a given particle, negligible rates of heat conduction between particles and within the fluid, constant properties, and negligible thermal expansion. Convective heat transfer between the phases is assumed to be proportional to the difference between the average phase temperatures at a given point. This model is mathematically given by ∂〈Tf 〉f hsf asf uD ∂〈Tf 〉f ∂t = ϵ4ρcp 5 <〈Ts 〉s - 〈Tf 〉f = - ϵ ∂z (3.12) f ∂〈Ts 〉s h a ∂t sf sf = - (1-ϵ)4ρc 5 (〈Ts 〉s - 〈Tf 〉f ) (3.13). p s Here, 〈Ts 〉s and 〈Tf 〉f are the volume-averaged solid and fluid temperatures, respectively, (𝜌cp) is the heat capacity (density and specific heat) of the given phase, hsf is the interstitial convective heat transfer coefficient, z is the position along the central axis of the bed, the subscripts s and f represent the solid and fluid phases respectively, and asf is the specific surface of the bed, defined as 87 Ap asf = Vp (1-ϵ) (3.14) [74] where Ap is the surface area and Vp the volume of an average particle in the bed. This model accounts only for the bulk movement of fluid and the convective heat transfer between the two phases. The second approach builds on the Schumann model by removing the assumption of negligible heat conduction within both phases. The solid phase is assumed to be completely continuous and axial heat conduction is accounted for with effective thermal conductivities. In the case of a packed bed, the contact resistance to conduction between particles is built into the effective thermal conductivity of the solid phase, allowing the non-agglomerated bed to be treated as a continuous phase. This is known as the Continuous-Solid-Phase (C-S) model [64], [75] and is given by ∂〈Tf 〉f keff,f ∂2 〈Tf 〉f hsf asf uD ∂〈Tf 〉f ∂t = ϵ4ρcp 5 ∂z2 + ϵ4ρcp 5 <〈Ts 〉s - 〈Tf 〉f = - ϵ ∂z (3.15) f f ∂〈Ts 〉s keff,s ∂2 〈Ts 〉s hsf asf ∂t = (1-ϵ)4ρcp 5 ∂z2 - (1-ϵ)4ρcp 5 (〈Ts 〉s - 〈Tf 〉f ) (3.16) s s Here, keff,s is the effective thermal conductivity of the solid phase over the cross-section of the REV domain (determined from effective conductivity models for packed beds), and keff,f is the effective thermal conductivity of the fluid phase over the cross-section of the REV domain, given as keff,f = ϵkf (3.17). The third approach removes the assumption of uniform temperature within a particle, choosing to treat the solid phase as a series of single spherical particles with concentric temperature profiles. The effect of axial fluid dispersion from flow through the tortuous paths within the pore space is also included via the axial dispersion coefficient Ddzz , and the convective heat transfer 88 between the phases is redefined to occur at the surface of each spherical particle. This model is known as the Dispersion-Concentric (D-C) model [64], [76] and is mathematically given by ∂〈Tf 〉f 1 k ∂2 〈Tf 〉f hsf asf uD ∂〈Tf 〉f ∂t = ϵ 14ρceff5 + Ddzz 2 ∂z2 + ϵ4ρcp 5 r2 ∂r ? (3.19) s ∂Ts -ks ∂r = hsf (Tsf - 〈Tf 〉f ) (3.20) where r is the radial position within the spherical particle and Tsf is the temperature at the surface of the spherical particle. For the purposes of this dissertation, the C-S approach was chosen because it includes axial conduction along the bed and is less complex than the D-C approach, as it avoids the need for measuring axial dispersion coefficients. 3.5. Summary This chapter summarized the relevant porous medium heat and mass transfer theory for single-phase flow through homogeneous porous media. The volume averaging treatment was qualitatively summarized, with the important mathematical laws presented and discussed. The fluid mechanics in isotropic, disordered media were discussed, highlighting the pore-scale Reynolds regimes and the lack of concern for macroscopic turbulent effects. Conduction heat transfer in packed beds was briefly discussed and a model for the packed beds of spheres was selected from the literature. The origins of the non-local-thermal-equilibrium equations for porous media heat transfer were highlighted. The three simplified 1D, single-phase flow, non- local-thermal-equilibrium energy models were presented and the Continuous-Solid-Phase modeling approach was selected. The material in this chapter forms the fundamental basis for experimental measurement of the interstitial convective heat transfer coefficients in packed beds of magnesium-manganese oxide. 89 CHAPTER 4: INTERSTITIAL CONVECTIVE HEAT TRANSFER IN PACKED BEDS OF MAGNESIUM-MANGANESE OXIDE This chapter overviews the experimental measurement of the interstitial convective heat transfer coefficients in packed beds of magnesium-manganese oxide and presents convection correlations that can be used for modeling the convective heat transfer in a reactive packed bed of magnesium-manganese oxide particles. 4.1. Review of the Experimental Literature The literature on experimental measurements of forced convection in packed particle beds is quite expansive and spans back many decades. A variety of both steady state and transient approaches have been taken. The steady state techniques involve some sort of heat generation or heat sink within the particle(s) of interest within a larger bed. The temperatures of these particles are measured locally, either at the particle surface or inside the particle. Some examples include (a) electrical heating of a single sphere in a bed of inert spheres and measuring the surface temperatures [77]-[79], (b) inductive heating of an entire particle bed [80], (c) constant drying of water-soaked spheres for simultaneous heat and mass transfer measurements [81]-[83], and (d) resistive heating of the particle bed with a direct current [84]. The transient techniques are different variations of responses to a step-change or an oscillatory frequency cycle in the inlet temperature [85]-[88]. A wide collection of studies is compiled by Barker [89] and is a useful reference to quickly compare much of the literature. His work compiles the surveyed literature into a table and series of charts, but does not attempt to generate a combined correlation. However, this is done in at least four instances in the works of Wakao et al [76], Gupta et al [90], Whitaker [Error! Reference source not found.] , and Sen Gupta & Thodos [91]. In the case of Wakao et al, the data was re-processed and modified to fit their specific 1D packed bed heat transfer model. The remainder simply generated 90 correlations of the data they surveyed, as published. Each general correlation is built from a different collection of studies, with some overlap between the four. In addition, there are a few theoretical studies that generated Nusselt correlations for packed beds [92],[93]. The Gnielinski correlation [93] in particular is a widely used theoretical relationship. Built on the equations for convective heat transfer from a flat plate, Gnielinski combines the asymptotic solutions for laminar and turbulent flow and, using the mean pore velocity and defining the length scale as the particle diameter, presents an expression for the heat transfer from a single sphere. This is extended to packed beds by multiplying the single-sphere equation by an empirical geometry factor [93]. Finally, there are also more recent studies that use direct numerical simulation of convective heat transfer within a packed bed of spheres to develop a Nusselt correlation, such as the work by Singhal et al [94]. A sample set of correlations from the literature have been collected and are compiled in a particle-diameter based Nusselt vs Reynolds number plot for comparison with experimental data in a later section. Further, the data collected by Gupta et al [90] have been scanned and compiled into a Chilton-Colburn factor [95] plot as an alternative means of presenting data. In this work, both a steady state and a transient approach are considered for measurement of the interstitial heat transfer coefficient (HTC) in packed beds. The remainder of this chapter develops the methods and presents the experimental results. 4.2. Steady State Measurement In general, the steady state approach proceeds as follows. First, a heat source or sink of some kind is applied to the particles in the packed bed. Next, a fluid flow is forced through the bed at a constant flow rate. The system is then left alone until reaching equilibrium when both phases exhibit stable, measurable temperature profiles across the bed. How the results are 91 extracted from these data depends on the modeling approach. Most measurements are structured such that a phase temperature difference can be directly computed, while other experimentalists have chosen to use the two-equation Continuous-Solid phase (C-S) model to extract interstitial HTCs. In these instances, the equations are derived in two dimensions, the heat source is applied along one of the boundaries of the porous material, and the HTC is extracted through an iterative modeling approach by finding the best-fit between the data and the simulations. The work by Garrity et al [96] provides a good example of this technique as applied to porous metal foams. The big advantage of steady state measurements is that corrections for heat loss are easily made. During an experiment, the porous material experiences a steady heat loss to the environment of some magnitude. Because it is steady, it can be compensated for in the data extraction process through calibration of the experimental structure or some other appropriate method. Additionally, steady state measurements eliminate the time-dependent term from the energy balance, removing the density and specific heat of the solid phase from consideration. This reduces the number of other measurements built into the process and could allow for studying interstitial convection within a porous material of unknown thermal mass. The big disadvantage is that distinguishing between the two temperature fields (solid and fluid phases) can be difficult to do with enough precision for an accurate interstitial HTC measurement. The porous material, particularly in the case of a packed bed, can be very dense with small pores. Ensuring that one is measuring the temperature of the expected phase requires precise placement of the temperature sensors within the bed and likely requires extra design efforts to isolate the sensors from the other phase. Depending on the geometry of the system (total volume, pore size, particle shape), it may require a large quantity of temperature measurements to fully resolve the two temperature fields. Insertion of probes into the bed at 92 multiple locations may also disrupt the fluid velocity profile downstream of each probe enough to alter the measurement. Finally, in some instances it may not be possible to distinguish between the thermal fields of the two phases: the thermophysical properties of the fluid-solid combination may cause thermal equilibrium to be reached so rapidly after the fluid enters the bed that there is not enough space to reliably capture the temperature profiles before equilibrium. In these cases, the steady state approach may be impossible. In fact, Wakao et al [76] go as far as to say that this approach is completely unfeasible for packed beds under any circumstances without the particles within the bed acting as a source/sink, generating or absorbing heat at steady conditions. In an attempt to avoid these problems, the direct current (DC) Joule heating approach is chosen, similar to that used by Glaser and Thodos [84]. This approach forgoes the C-S model & iterative optimization approach and extracts the interstitial HTC using only a few energy balances and direct temperature measurements. However, achieving accurate results with this approach with a reasonable confidence level proved to be very difficult and was subsequently abandoned for further consideration in this work. Thus, the work done towards steady state measurements is discussed separately in Appendix C. 4.3. Transient Measurement The basics of the transient approach are as follows. The packed bed of particles is taken from some initial thermal equilibrium to some final thermal equilibrium by a fluid flow, the transient response of the system is recorded, and the resulting data is used to extract the interstitial HTC. One specific approach to transient measurements is the single-blow technique. A type of step-response technique, it was first proposed for testing heat exchangers [97]-[98] and has since evolved to be applicable to various porous media [99]-[104]. Particles are packed into an 93 enclosure of some geometry (Cartesian or cylindrical) and supported on either end with some structure such that the bed is well-packed and remains fixed, yet fluid can flow through the bed unobstructed. The bed is brought to some initial starting state, either a “hot” or “cold” state depending on how the measurement is to be performed. Ideally, this is a uniform bed temperature. Once equilibrium is reached, a fluid of known temperature and flow rate is forced through the bed, bringing it to some final equilibrium state over a length of time. Temperatures of the fluid at the inlet and outlet and the pressure drop across the bed are recorded until thermal equilibrium is reached. The system is then modeled using one of the two-phase non-equilibrium porous media continuum models with the transient experimental data as inputs. The outlet temperature of the experimental system is modeled for repeated iterations of the interstitial HTC until a best-fit is reached between the modeled and measured outlet temperatures. The resulting HTC is considered to be the correct value. The biggest advantage of the transient approach is that experimental data are relatively easy to acquire, provided the thermophysical properties of both the fluid and solid in the packed bed are known. The fluid temperature across the face of the inlet and the outlet, the fluid flow rate, and the pressure drop across the bed are the only required datasets. These measurements are relatively simple to make with reliable instruments. However, this approach becomes more complicated if the thermophysical properties of the packed bed are uncertain. The equations require that the density, specific heat, and effective thermal conductivity of the packed bed be known. Additionally, all energy flows into and out of the test section structure must be accounted for (or be negligible) to get accurate results. This can create problems if the scale of the test section structure is small enough that the transient heat transfer between the packed bed and the wall during the experiment observably impacts the outlet fluid temperature data. Further 94 difficulties can arise when trying to reduce the dimensionality of the problem. A 1D analysis is ideal but requires a Darcian velocity profile, uniform inlet and outlet temperatures, limited channeling of the fluid through the bed, and no significant radial temperature gradients within the packed bed. All of these factors must be taken into account in designing the experimental setup. 4.3.1. Experimental Setup and Method A custom setup was built for performing the transient measurements. Three stock 25% glass-filled PTFE tubes from McMaster-Carr were machined to make a test section and two flow-development sections that could connect to standard pipes and tubing, support a vertically- oriented packed bed, and minimize the thermal mass of the tube wall along the test section. 25% glass-filled PTFE was chosen for its machinability, low thermal mass in comparison to metal tubing, suitably high melting temperature (260 °C), low thermal conductivity (0.45 Wm-1K-1), suitable rigidity, and ease of acquisition at reasonable prices compared to other plastics such as PEEK. To make the test section, one stock tube was turned to an external diameter of approx. 7.6 cm (3”). The interior diameter (nominally 2”) was bored on the lathe with a nominal 2” diameter drill bit to straighten out the internal diameter. The final diameter was measured with a Vernier calipers at 50.6 mm (1.99”). The thickness of the tube wall along a length of 15.2 cm (6”) was reduced to the thinnest extent possible until further machining risked puncturing the wall at any locally-thin spots or ripping the part entirely due to torsion at the base. The final wall thickness ended up at approx. 1.9mm. This thickness was achieved by placing a thick-walled aluminum tube inside the PTFE tube to provide support during the final turning and smoothing process. The test section was then cut to length and faced on either end to create flanges. Six ¼” holes were drilled for bolts in each flange at even intervals. Additional holes were drilled through the 95 exterior flange walls in the transverse direction, tapped, and fitted with stainless steel compression fittings to allow access for pressure sensors, thermocouples, and a vent at the inlet. The two additional PTFE tubes were machined to an OD that matched the OD of the test section flanges. These PTFE tubing parts form the core of the experimental setup. The setup is assembled as follows. The upstream flow-development section is packed with 6mm plastic spheres, while the downstream is packed with 6mm porcelain spheres to avoid any melting concerns from continual exposure to hot air. Each flow-development packed bed is permanently fixed using ¼” thick Nomex honeycomb disks held in place by screws drilled into the tube wall. These parts are independent of the test section and are bolted to either flange. This is done to (a) homogenize the velocity profile such that plug flow is achieved, (b) minimize the entrance and exit effects of the test section, and (c) assist with developing uniform temperature profiles across the tube face on either side of the test section. Aluminum mesh supports installed in the test section support the packed bed. The outlet mesh is fixed to the inside wall of the outlet flange using short self-tapping screws. The inlet mesh is a removeable disk that is wedged in place after packing the bed. Three K-type thermocouple probes from Omega Engineering are installed at different radial locations just downstream of the support mesh at the outlet to record the gas temperature. Preliminary tests indicated that the temperature profile at the inlet was uniform, so only a single probe was placed just upstream of the inlet along the centerline. This was beneficial because it made it easier to change the packed bed and insert/remove the upstream support mesh. The probes were calibrated prior to installation to a Pt-100 high-precision thermometer within a hot water bath between 22- 75 °C. Additional points were taken using ice water and boiling water to expand the calibration range. The Inconel sheathing and mica insulating powder that covered the thermocouple 96 junctions were carefully removed from the probes using a Dremel to expose the bare wire. This was done to remove any transient delay errors caused by the sheathing acting as a thermal barrier between the air flow and the thermocouple junction. Two Dwyer Instruments Series-626 pressure transmitters record the air pressure at the inlet and outlet. High-temperature silicone rubber gaskets of medium durometer (50A) and thickness 0.64 cm (¼”) seal the system. Additional sealing around the tapped holes is done as- needed with Permatex Ultra-Black brand high-temperature silicone rubber motor gasket sealant. Any small leaks are identified with a soap-solution leak check at the highest pressure observed during all experiments and subsequently sealed with the silicone sealant. Alicat MC-Series mass flow controllers regulate the air flow rate. Additionally, the precise flow rate for each experiment is determined by performing a pre-check with the Alicat M-Series flow meter of range 0-500 SLPM. The check is done every time the system is reconfigured with a new packed bed to ensure accurate quantification of the flow rate. Finally, E-type thermopiles installed in the tube wall at several locations along the test section measure the temperature gradient across the wall during each experiment to quantify the transient heat transfer effects between the bed and the test section structure. An illustration of the test section made by Polo Design Services is shown in Figure 4.1. 97 Figure 4.1. Transient Experimental Setup, Test Section. Illustration of the test section in the transient interstitial heat transfer coefficient experimental setup. 98 The working fluid for all the experiments is compressed air from the main building air compressor that is fed through a pressure regulator set at 60 psi to remove any oscillations in the flow rate due to tank pressure variations, then fed through a combination oil trap/particle filter/water filter setup to dry and purify it prior to entering the mass flow controllers. A set of flexible ¾” tubing and valves are placed throughout the system to allow for the upstream section to be kept cool during the pre-heating step, for the pre-heating to be done with a counter flow of hot air, and for the blow-down step to be done using a top-to-bottom flow of cool air. A 1/4” NPT valve is installed even with the face of the packed bed in the inlet flange for venting hot air after passing through the bed. Air is heated by a 400 W 316 stainless steel in-line duct heater from Omega Engineering. Air temperature is controlled by connecting a JLD-612 PID temperature controller to the centerline outlet thermocouple. Temperature is set at 72 °C such that each blow-down experiment occurs over an approx. 50 °C range, given that the average inlet air temperature from the compressed air system is approx. 22 °C on average. An illustration of the flow management system from Polo Design Services is shown in Figure 4.2. 99 Figure 4.2. Transient Experimental Setup. Illustration of the entire transient interstitial heat transfer coefficient experimental setup. A small amount of fibrous high-temperature aluminosilicate insulation from Unitherm is placed inside the flange on the inlet end of the test section to reduce losses to the cooler upstream length so that the initial bed temperature is as close to uniform along the axis as possible. Radial heat loss is prevented by insulating the entire test section and flanges with very thick-walled (approx. 3.5” thick) fiberglass pipe insulation wrapped in a foil-lined paper jacket. Additional fibrous aluminosilicate insulation then fills any small gaps between the structural components and the insulation. The hot air pipeline is insulated in a similar fashion using appropriately sized insulation jackets. The entire setup is vertically oriented to prevent any settling of the particles and subsequent channeling of the air during the experiments. The setup is supported by a 1” square slotted framing structure. Support brackets are placed underneath the most-downstream PTFE flange where the PTFE parts connect to the aluminum tubing and are mounted to the frame. Additional braces are placed vertically along the downstream PTFE flow-development 100 part to prevent any buckling of the column. The top portion is braced laterally with two framing brackets, but allowed to expand vertically when heated. The experimental procedure is as follows. The bed is heated with a flow of 50 SLPM from the main flow controller, through the heater, upward through the bed, then vented through the ¼” valve at the inlet. 50 SLPM is chosen so that some level of turbulent mixing occurs within the hot air supply line to give a radially-uniform temperature profile so that the bed heats evenly while maintaining a low enough flow rate that the flow does not choke when exiting through the ¼” NPT valve at the inlet. Simultaneously, the secondary flow controller is set at 5 SLPM and allowed to flow through the upstream sections to the inlet, then out the ¼” valve along with the hot air. This keeps the upstream sections cool on the inside during the pre-heating in order to have a uniform inlet temperature profile. The system is allowed to heat until a steady-state is reached across the bed (approx. 1.5-2 hours depending on the bed material). At this point, the cool upstream flow is stopped. The main flow controller is vented to ambient by switching a series of valves. The setpoint is changed to the experimental flow rate and the device continues to vent to atmosphere for a short time while the controller stabilizes (approx. 5-15 sec depending on flow rate). During this time, more valves are switched to cut off the hot air line, the secondary flow controller, and the ¼” valve vent such that all air will flow through the test section and then vent to ambient. The heater is also turned off. Finally, the last valve is switched and the main flow is forced through the bed, creating the inlet temperature step-input to affect the blow-down transient response of the bed. Temperature data across the face of the inlet and outlet and pressure drop across the bed are recorded until the inlet and the outlet read approx. the same temperature. Data is acquired using a Labjack 24-bit T7-Pro data acquisition board at half- second increments. 101 Ten datasets are collected for a given packed bed at flow rates between 50-500 SLPM at increments of 50 SLPM each. Because the porous media theory assumes perfect spheres, packed beds of 3, 5, and 8mm 304SS and 4 & 6mm non-porous alumina spheres are tested to form the baseline correlation and validate that the setup behaves as expected. A precise estimate of the diameter and density of each set of spheres was found by taking a random sample of 20 particles, measuring the diameter and mass of each, then calculating the average diameter and average density. These values were in good agreement with those expected and are used as inputs for the model. Two sets of 1/1 Mn/Mg particles are also tested under the same conditions: (1) granular particles with a median effective diameter of 3.7mm and (2) cylindrical particles of diameter 3.5mm and height 2.3mm (effective spherical diameter of 2.98mm [74]). Average dimensions and densities of both particle sets were quantified as a combined effort between the author, Dr. Kelvin Randhir, Phillip Schimmels, and Michael Hayes. 4.3.2. Thermal Effect of the Tube Wall and Volume-Averaged Model It is important to make note of the scale chosen for these experiments. Constraints with the amount of 1/1 Mn/Mg material available, limitations of the lab compressed air supply, and both the cost and availability of the structural components dictated that a small-scale setup was more feasible than a larger scale one. However, the preliminary designs and subsequent modeling analysis made it clear that at this scale, the energy effects of the test section had to be accounted for to achieve any kind of reasonable results using this measurement technique. Figure 4.3 gives an example. 102 Figure 4.3. Transient Wall Heat Transfer Effect. Example plot showing the extra energy transferred from the wall to the system during an experimental run. The median experimental case of 5mm particles with a flow rate of 250 SLPM is used for this example. The measured average outlet temperature is shown along with three outlet temperature curves simulated using the C-S model for three different interstitial heat transfer coefficients (HTC). The three simulated outlet curves show that variation in only the interstitial HTC without any compensation for the transient wall heat transfer effect is incapable of closing the energy balance between data and simulation. The simulated outlet curves in Figure 4.3 come from the C-S model for three different interstitial HTCs. The figure shows that a single-parameter optimization using the C-S model without any additional terms is incapable of achieving a good fit with the data and thus is unsuited for use in its basic form. Further, because the experimental data take noticeably longer to cool than the three simulations, it is clear that the data contain more energy than the model is accounting for. Thus, the energy balance is not closed between the model runs and the experimental data, and a best-fit result cannot be confidently achieved. Because confidence in 103 the thermal mass of the bed and the magnitude of the air flow rate is high, the only remaining source of the extra energy is the surrounding experimental structure. Thus, a variation of the C-S porous medium heat transfer model is derived from an energy balance to include the effects of the wall in the analysis. The system of equations used for the data analysis is as follows. The solid phase energy equation is given by, ∂〈Ts 〉s ∂2 〈Ts 〉s (1-ϵ)ρs cps = keff,s - hV <〈Ts 〉s -〈Tf 〉f = (4.1) ∂t ∂z2 with boundary conditions ∂〈Ts 〉s ∂〈Ts 〉s -keff,s ∂z @ = keff,s ∂z @ =0 (4.2). z=0 z = Lbed Here, keff,s is the effective thermal conductivity of the packed bed, hV is a volumetric heat transfer coefficient given by hV = hsf asf (4.3) where hsf is the interstitial heat transfer coefficient and asf is the specific surface of the packed bed. 𝜖 is the bulk porosity of the packed bed, i.e., the phase volume fraction of the fluid phase. Because there are only two phases, the volume fractions of each phase must sum to 1. Thus, bulk porosity can be defined using the definition of the volume fraction given by Eq. 3.1 and is determined by the following equations: Vs + Vf = VCV (4.4) Vf ϵf = VCV (4.5) Vs Vf ϵs = VCV =1- VCV = 1 - ϵf (4.6) ϵf = ϵ, ϵs = 1 - ϵ (4.7). 104 Bulk porosity is calculated from Eq. 4.4-4.7 using the known internal volume of the test section containing the packed bed and the volume of solid particles (determined from the mass of the particles and the density of the solid phase). The fluid phase energy equation is given by, f ∂〈Tf 〉 ∂2 〈Tf 〉 f ∂〈Tf 〉 f 2 ϵρf cpf ∂t = ϵkf ∂z2 + hV <〈Ts 〉 s - 〈Tf 〉 f = - m"f cpf ∂z +R q"w (4.8) bed where m"f is the mass flux of air through the bed (calculated from the fixed air flow rate and the cross-sectional area of the test section), Rbed is the internal radius of the test section, and q"w is the heat flux between the wall and the air and is calculated from the measured temperature gradient across the tube wall. The assumption is that the contact points between the particles and the tube wall are small enough that all wall energy is transferred to the fluid phase. q"w is calculated from the measured temperature gradient data using a Fourier’s law approximation as kw ΔTw q"w = - tw (4.9) where kw is the thermal conductivity of the wall (25% glass-filled PTFE), tw is the average tube wall thickness (approx. 1.9mm), and ΔTw is the measured temperature gradient across the tube wall. The tube wall is thin enough that tw is assumed to be sufficient instead of an ln(r2/r1) length scale term. The boundary conditions are given by ∂〈Tf 〉f 〈Tf 〉f (z = 0) = Tin , -kf @ =0 (4.10). ∂z z = Lbed This modified C-S model will be used to extract the interstitial HTCs from the transient experimental data. The derivation of the model can be found in Appendix D. 4.3.3. Data Analysis: Numerical Method, Data Extraction, and Uncertainty To extract the interstitial HTCs, the modified C-S model is solved using a 1D explicit finite difference numerical scheme. Diffusion terms are discretized with a second-order accurate 105 central difference scheme and advective terms with a second-order accurate upwind scheme. Time integration is done with a two-step Adams-Bashford approach. The spatial domain is divided into 31 nodes and a small time-step of 1e-4 s is used for all experiments except those at 50 SLPM, when a time-step of 5e-3 s is used. The inlet temperature boundary condition is prescribed as a function of time by the inlet temperature data. The remaining Neumann boundary conditions are discretized using the ghost-node approach assuming zero flux across the boundaries. The initial condition is assumed to be a uniform temperature for both phases across the domain and is taken to be the average of the recorded inlet and outlet temperatures after the bed reaches steady state at the end of the pre-heating step and prior to beginning the blow-down step. Thermophysical properties of the air are allowed to vary with temperature and are calculated from curve-fits of the data at 1 atm between 0-100 °C [105]-[107]. Density of air is calculated at each node and time step assuming the ideal gas law. The specific heat and effective thermal conductivity of the solid phase are calculated at the midpoint temperature of the blow- down step (approx. 48 °C). Specific heats of the 304 stainless steel and alumina particles are well-established in the literature and are used in the model. Effective thermal conductivities of the packed beds of spheres are estimated using the Hadley General Average model. For the 1/1 Mn/Mg particle beds, the low-temperature specific heat measurements in Appendix E and the experimentally measured effective thermal conductivities of each of the particle beds [30] are used. The bulk porosity of each particle bed is calculated from the mass of the bed, the density of the solid material, and the volume of the cylindrical test section. The inlet temperature, inlet pressure, outlet pressure, and wall temperature gradient data for each experimental run are curve-fit using a set of piece-wise 6th order polynomials. A pre- 106 processing routine prepares the data for input into the model such that a value can be calculated for each discrete point in space and time during the simulations. The points in time defining the piece-wise split for each function are identified prior to simulations and checked to give good representations of each dataset in Excel, then input as constants into the pre-processing routine. Curve-fitting is done in MATLAB using the polyfit functions [108]. Pressure drop across the bed is assumed to be linear and is calculated at each node from the curve-fits of the inlet and outlet pressure signals. Flow conditions are defined as a constant mass flux of air through the bed, calculated from the recorded flow rate. Data is partitioned prior to simulations such that time t = 0 is the point in time at which the pressure signals first indicate the presence of the cooling air flow, ensuring that no data during the pre-heating step is considered. The best-fit volumetric HTC is determined by an iterative single-parameter optimization using the fminsearch function from the MATLAB Optimization Toolbox [108]. The function minimizes the normalized root-mean-squared error (NRMSE) between the simulated and measured outlet temperatures, calculated by N 2 9∑i=1)Texp(i)-Tmod (i)* N NRMSE= mean(Texp ) (4.11). Tolerances of 5e-4 for the NRMSE and of 500 W/m3K for the volumetric form of the HTC used in the model were used. The initial guess for each optimization was calculated from the Wakao et al correlation [76] and the specific surface asf for a given packing. The resulting optimized hV is then converted to the surface-area dependent hsf by Eq. 4.3. A secondary consequence of the transient wall heat transfer effect is that the inlet and outlet temperature measurements reach a long-time quasi-equilibrium state where the slow removal of the remaining heat from the test section causes a small difference between the inlet 107 and outlet such that the outlet signal has a very long tail. This posed an obstacle to data extraction. To do the optimization analysis, some cutoff point in time had to be defined consistently across all the datasets beyond which no more data was considered by the optimizer. The natural choice is the point at which the difference between the average inlet and outlet temperatures drop within the maximum uncertainty range of the thermocouples. However, because of the long-tail effect, this point is not reached for a very long time and consideration of the tail introduces a heavy skew towards the tail in the optimization algorithm and the resulting HTC. Thus, a cutoff point is defined using the following criteria. A dimensionless time parameter is defined based on the thermal mass of the bed and the accumulated thermal mass of air, given as t ∫0 ∞ ṁ air cpair dt t = * 4mcp 5 (4.12) bed where t∞ is the time at which the experiment is considered complete, ṁ air is the mass flow rate of air used in the given experiment, andk p ? (4.22). air 4.3.4. Results The results for the transient interstitial HTC measurements are provided below. Figure 4.4 contains plots of (a) the dimensionless datasets and (b) the dimensionless “best-fit” outlet temperature simulations that correspond to the measured interstitial HTC for the packed bed of 3mm 304SS spheres. The dimensionless time used is the same used for determining the end- point of each experiment, as described by Eq. 4.12. Dimensionless outlet temperature is defined as Tout - T∞ θout = T0 - T∞ (4.23) where T∞ is the end-point equilibrium (i.e., the constant inlet temperature), and T0 is the initial temperature of the system prior to the blow-down step. Figure 4.4 and similar plots for the other six particle beds tested are included together in Appendix F. 111 Figure 4.4. Data and Simulated Best-Fit Curves. Dimensionless plots of the (a) measured and (b) best-fit simulated outlet temperatures for the packed bed of 3mm 304SS spheres. All ten runs (50-500 SLPM) are contained in each plot. The banded nature comes from longer time exposure to the transient wall heat transfer effect and skews lower flow rates towards larger error in the resulting interstitial HTC measurement than at higher flows. (a) (b) 112 In general, all ten datasets between 50-500 SLPM for a given particle size align reasonably well, but exhibit a banded structure. The general trend is that the curves to the right side of the band are at lower flows, then as the flow rate increases, the curves shift to the left of the band and become more closely aligned. This occurs in both the measured data and the “best- fit” simulated outlet temperature that corresponds to the optimal interstitial HTC for a given run. This banded nature comes from the transient wall heat transfer effect and is due to the amount of time each experiment takes to run. As the flow rate decreases, both the residence time of the air in the test section and the total time to complete the experiment increase. Because of this, lower flow experiments allow more time for heat absorbed by the structure during the pre- heat to diffuse into the test section and propagate into the data through convection between the air and the test section wall along the bed. Thus, when made dimensionless, the measured and simulated outlet temperatures do not line up perfectly because there is more energy removed overall in the data at lower flows than at higher flows. However, the quality of fit improves as the flow rate increases and the relative noise of the wall heat transfer effect is less. This is seen in Figure 4.5, where plots of three of the curves shown in Figure 4.4 for the bed of 3mm 304SS spheres at flow rates of 50, 250, and 500 SLPM are shown individually. 113 Figure 4.5. Increased Error at Lower Flow Rates. Plots of the measured and simulated “best-fit” outlet temperatures for the 3mm sphere bed at (a) 50, (b) 250, and (c) 500 SLPM. The observed overshoot & temperature-spike effect seen in the simulated data, an artifact of the error introduced to the simulation from the point-measurements of the transient wall heat transfer effect, becomes less pronounced as the flow rate increases. (a) 1.2 1 0.8 Data Model θout (-) 0.6 0.4 0.2 0 0 1 2 3 t* (-) (b) 1.2 1 0.8 Data θout (-) 0.6 Model 0.4 0.2 0 0 1 2 3 t* (-) 114 Figure 4.5 (cont’d). (c) 1.2 1 0.8 Data Model θout (-) 0.6 0.4 0.2 0 0 1 2 3 t* (-) The quality of the fit with the data clearly improves as the flow rate increases and the noise of the wall heat transfer becomes less significant relative to the volumetric heat transfer signal. In particular, the plot of 50 SLPM is noteworthy because of the spike in the simulated outlet temperature. This occurs in the simulation because of the incorporation of the wall temperature gradient (and associated heat flux) measurements. Removal of those measurements completely eliminates any spike-like behavior in the simulations. At lower flows, it seems that the thermal interactions (losses and gains) between the test section walls and the packed bed are not captured as well with the 3-4 point-measurements of the wall temperature gradient than at higher flows. This could mean that a more accurate measurement requires more points of wall temperature gradient measurements, or it could just be a product of the higher percentage of experimental uncertainty that these measurements introduce at lower flows relative to that at higher flows. Regardless of the precise cause, the consequence is that even after incorporating measurements of the transient wall heat transfer effect into the interstitial HTC measurement, 115 there is less confidence in the lower flow rate measurements than in the higher flow rate measurements. Figure 4.6 plots the optimal results as a Nusselt vs Reynolds number plot using the definitions given by Eq. 4.16 and 4.17. Figure 4.6. Results, Nusselt. Results from the transient experiments as Nusselt numbers covering the full experimental range with Reynolds numbers of 50-2000. The results for the seven different particle beds examined are shown along with plots of the compiled literature correlations for comparison. 180 Wakao/Kagui/Funazkri (1979) 1/1 Mn/Mg Cylindrical, 3.5x2.3mm 1/1 Mn/Mg Granular 3.7mm 160 Alumina, 6mm 304SS, 8mm 304SS, 5mm 304SS, 3mm 140 Alumina, 4mm Gunn (1978) Furnas (1930) Glaser & Thodos (1958) 120 Singhal et al (2016) (DNS Spheres) Achenbach (1995) Gnielinski (1978) 100 Gupta, Chaube, Upadhyay (1974) Whitaker (1972) Nudp(-) Sen Gupta & Thodos (1962) De Acetis & Thodos (1960) 80 Malling & Thodos (1966) McConnachie & Thodos (1963) 60 40 20 0 0 500 1000 1500 2000 Redp (-) From the figure, the following results are observed: 1. The measurements for the 3, 4, and 5mm diameter spheres align very well. These results match despite the difference in materials (304SS and alumina), meaning that they are independent of the thermophysical properties of the solid bed, as expected. 116 2. The measurements for the 6 and 8mm spheres have slightly smaller magnitudes than the 3/4/5mm, with that of the 8mm being slightly less than the 6mm. 3. The error bars for all measurements maintain approximately the same lengths as the flow rate increases. However, this corresponds to the error at low flows being a larger percentage of the measurement than at higher flows. This confirms the notion that the lower quality of fit at lower flow rates corresponds to a larger percentage of error introduced by the transient wall heat transfer effect. 4. The datasets for both the median 3.7mm diameter granular and the 3.5 x 2.3mm cylindrical (effective particle diameter of 2.98mm given by Eq. 4.15) follow distinctly different trends than the perfect spheres of 304SS and alumina. This was expected for the cylinders, as there is evidence in the literature that Nusselt numbers for packed beds of cylinders are lower than for beds of spheres because of a difference in packing arrangements [89]-[91]. However, the reasoning is less clear for the granular particles. These particles are granule in shape, i.e., they do not have a perfectly-defined or uniform geometry, but rather are variously ellipsoidal or spherical with very rough and unique surface variations. Further, because of their manufacturing process, they do not have a uniform diameter but rather a distribution of effective spherical diameter. These geometric effects mean that (a) the packing and porosity will be less homogenous at a local level than for perfect spheres and (b) it’s possible that locally more or less dense packing around one of the thermopiles across the wall could cause a change in the measured wall temperature gradient that causes the approximation to reflect the true effect less well than it may have for a more uniform packing. Regardless, the result is that these data do not line up well with the perfect spheres and should not be included in a 117 general correlation, but rather be given specific correlations that are only relevant to the specific material, particle shape, and diameter distribution. Further investigation of this behavior would be warranted in a future study to determine the precise causes. It may also be worth repeating these experiments in a different setup where the wall heat transfer effect is less prevalent to see if the results hold. It’s clear from the results that there is the possibility of a general correlation for the 3/4/5/6/8mm spheres. However, the slight decrease in magnitude of the Nusselt number for the 6 and 8mm spheres must be rationalized. The reason lies in the ratio of tube-to-particle diameter (D/dp). When packed in a tube, there is a region near the tube wall that has an oscillatory porosity as the packing becomes denser further into the core of the bed [66]. This occurs in every system and is lumped into the average bulk porosity calculation. However, there is evidence in the literature that when the D/dp ratio goes below 10 for a given system, the less-dense packing near the wall regions becomes substantial enough to cause non-uniformity in the flow profile and subsequently begin to break down the 1D plug flow assumption made by the porous media heat transfer models [80],[11010]. These effects get lumped into the interstitial HTC during experimental analysis. Table 4.2 shows that for the 6 and 8mm spheres, the D/dp ratio is less than 10, while for the rest of the spheres it is greater than 10. The average bulk porosity also increases by a small amount for the 6 and 8mm spheres relative to the 0.39-0.4 exhibited by the 3/4/5 mm spheres. 118 Table 4.2. Tube-to-Particle Diameter Ratios. Effective particle diameter (dp), tube-to-particle diameter ratio (D/dp ), and bulk porosity for each packed bed of spheres. D/dp becomes less than 10 and the bulk porosity increases for the 6 and 8mm diameter sphere beds. Material dp (mm) D/dp (-) 𝜖 (-) 304SS 7.99 6.3 0.43 304SS 5.00 10.1 0.39 304SS 3.03 16.7 0.39 Alumina 4.20 12.0 0.40 Alumina 6.05 8.4 0.42 If the geometry differences are accounted for in the correlation, it may be possible to collapse the five datasets to a single relationship. When the data are plotted as a Chilton-Colburn factor with bulk porosity included in the y-axis, the five datasets do indeed collapse. Figure 4.7 shows all the transient data plotted in this fashion. The Chilton-Colburn factor is defined as 1 -1 - jH =Nu Re Pr 3 (4.24) where Pr is the Prandtl number for air and the Reynolds and Nusselt numbers are given by Eq. 4.16 and 4.17. 119 Figure 4.7. Results, Chilton-Colburn. All interstitial HTC data plotted as a Chilton-Colburn jH factor with porosity variation incorporated into the y-axis. Addition of the porosity factor collapses the data for all five beds of perfect spheres. The 1/1 Mn/Mg particles continue to follow unique trends. 0.1 !jH (-) Scanned Data, Gupta et al (1974) 1/1 Mn/Mg Cylinders, 3.5x2.3mm 1/1 Mn/Mg Granular, 3.7 mm Alumina, 6mm Alumina, 4mm 304SS, 3mm 304SS,8mm 304SS, 5mm 0.01 50 500 5000 Redp (-) The figure shows that the distinctly different trends observed in the Nusselt plot continue for the 1/1 Mn/Mg particles. The cylinders have become more aligned with the spheres, but still lower in magnitude, and the granular particles continue to follow a unique trend. Further, the measurements at 50 SLPM for both sets of 1/1 Mn/Mg particles deviate from the general trend of the literature data surveyed and correlated by Gupta et al [90], indicating that these particular measurements could be unreliable. However, the incorporation of the porosity differences into the trend do indeed collapse the datasets for the perfect spheres sufficiently to generate a general correlation for the interstitial convection in packed beds of spheres. However, the skew in datapoints and size of the error bars at lower flow rates shows that this relationship will have a larger error at Reynolds numbers below 250, and certainly should not be used at Reynolds 120 numbers below around 75 without the addition of more accurate low Reynolds number datapoints to the correlated data. The correlation for packed beds of uniform spheres is given by Eq. 4.25 as ϵjH = 0.2733 Re-0.3194 dp (4.25). This relationship was fit using the MATLAB Curve-Fitting Toolbox with an R2 value of 0.955 and is valid for packed beds of spheres for particle-diameter based Reynolds numbers between approx. 75-2000. Figure 4.8 shows Eq. 4.25 and the data for beds of spheres for visually comparing the fit with the data. Figure 4.8 Results, Chilton-Colburn, Spheres. Chilton-Colburn factor plot of just the data for the perfect spheres, with Eq. 4.25 shown for comparison. Errors become less significant as the Reynolds number increases. The larger errors in the low Reynolds number measurements come from a larger error introduced by the transient wall heat transfer effect at lower flow rates. 0.1 Alumina, 6mm Alumina, 4mm 304SS, 3mm 304SS,8mm 304SS, 5mm Power Law Correlation !jH (-) 0.01 50 500 5000 Redp (-) 121 Finally, even though the data for the 1/1 Mn/Mg particle beds do not collapse into a general correlation, it is still useful to generate Nusselt correlations that can be used for systems that use these specific particles. The correlation for the 3.7mm granular particles is given by, Nudp = 0.2361 Re0.8557 dp (4.26) with an R2 value of 0.989. The correlation for the 3.5 x 2.3mm cylindrical particles is given by, Nudp = 0.2869 Re0.7598 dp (4.27) with an R2 value of 0.987. Both equations are plotted with the datasets as functions of Reynolds number in Figure 4.9. Figure 4.9. Results, Nusselt, 1/1 Mn/Mg Particles. Nusselt vs Reynolds plot of the data for the 1/1 Mn/Mg particles, with the correlations shown for comparison. The power laws slightly overpredict the heat transfer at low Reynolds numbers and would benefit from acquiring more accurate low Reynolds number measurements. 122 The figure shows that these relationships will likely mildly overpredict the convective heat transfer at Reynold numbers below 75 and should either be restricted to the ranges of data (50- 1000 for the granule and 50 – 800 for the cylindrical) or should be used in low Reynolds number regime systems only when the system has a low sensitivity to changes in the interstitial convection heat transfer coefficient. 4.3.5. Lumped Thermal Analysis of Particles Use of the modified C-S model for extracting the interstitial HTC results included the assumption that there were no significant temperature gradients within the individual particles in the packed bed. To validate this assumption, a lumped thermal analysis is done for each of the seven particle types tested. A Biot number for a single sphere undergoing convection heat transfer is defined as hsf rp Bi = kp 3 (4.28) where hsf is the measured interstitial HTC for a given particle type and Reynolds number, rp is the radius of the particle defined as half of the effective particle diameter dp, and kp is the thermal conductivity of the pure solid material. Literature values were found for the thermal conductivity of 304 stainless steel and alumina in the purity range of the particles used. The effective thermal conductivity model by Hayes et al [30] for the 1/1 Mn/Mg particle beds includes a calculation for the thermal conductivity of the 1/1 Mn/Mg material and is thus used for this analysis. The resulting Biot numbers for each measured datapoint are shown in Figure 4.10. 123 Figure 4.10. Lumped Thermal Analysis. Lumped thermal analysis for an individual particle in the packed bed. The Biot numbers for the spheres are very low and thus intra-particle temperature gradients during the interstitial HTC analysis are of no concern. Those for the 1/1 Mn/Mg particle beds are somewhat larger, trending above the accepted threshold for lumped analysis of Bi = 0.1, yet are still relatively small. This is attributable to a low thermal conductivity of the solid material at the median experimental temperature of 48°C (approx. 0.94 W m-1 K-1 [30]) . It is possible that small temperature gradients within the 1/1 Mn/Mg particles could cause a portion of the observed interstitial HTC deviations from the general trend for the spheres. The figure shows that the Biot numbers for all five particles of 304SS and alumina spheres are very small, well below the accepted threshold of Bi = 0.1 for the validity of a lumped thermal analysis. Thus, for these particles, intra-particle temperature gradients are of no concern and the assumption of uniform particle temperature made by the modified C-S model is valid. However, for the two 1/1 Mn/Mg particles, the Biot numbers are as large as 0.4 for the granular and 0.2 for the cylindrical particles. The thermal conductivity of the 1/1 Mn/Mg material at the 124 median experimental temperature of 48°C calculated from the model of Hayes et al [30] is much smaller than that of the other two materials (approx. 0.94 W m-2 K-1 vs 14.4 and 18 W m-2 K-1 for 304SS and alumina, respectively), which becomes evident in the different Biot number trends in Figure 4.10. These values are above the 0.1 threshold, but are also not very far above it. Therefore, it is possible that some portion of the deviation observed in the interstitial HTC results for the 1/1 Mn/Mg particles from the trend of the beds of spheres could be explained by small temperature gradients within the individual particles. However, because the Biot numbers are still quite small, the geometric effects are likely to have a much larger impact on the observed deviations from the general trend of the spheres. 4.3.6. Scale Analysis of the Transient Wall Heat Transfer Effect By far the most challenging element of designing a setup for transient measurement of the interstitial HTC in packed beds was properly handling the transient wall heat transfer effect during the experiments. Great lengths were taken to resolve this problem, and it is the author’s belief that there is not much room for improvement on the approach chosen to measure the effect in this work. As mentioned before, the scale of the test section was originally chosen over concerns of the availability of enough material for the packed bed, the capacity of the available compressed air system, and the cost of flow control devices for regulating the very high flow rates that a larger diameter setup would need in order to reach higher Reynolds numbers. However, the difficulty of compensating for and the error introduced by the transient wall heat transfer effect were substantial enough that a different set of design decisions likely should have been taken. It is thus of great interest to know under what conditions will this effect vanish from the data such that no transient wall heat transfer measurements and modifications to the porous media models are required. 125 The core problem is that the signal-to-noise ratio of the system is too small, if the heat transfer from the packed bed is considered as the “signal” and the heat transfer from the tube wall considered as the “noise”. Physically, the signal depends on the inside volume of the tube containing the test section, while the noise depends on the total area of the inside surface. Thus, it follows that the noise will vanish as the ratio of bed radius to bed length (R/L) increases. The question is thus “For a given Reynolds number, at what R/L ratio does the transient wall heat transfer become negligible?” To this end, a parametric study is done with a dimensionless form of the modified C-S model used in this work. The scaling parameters and dimensionless equations are as follows. Time is scaled by a parameter similar to a bed Fourier number, given by keff,s t* = t (4.29) (1-ϵ)ρs cp L2bed s where Lbed is the length of the packed bed along its axis. Space is scaled by taking the bed length as the characteristic length scale and defining a dimensionless length as z z* = L (4.30). bed The solid and fluid temperatures are scaled by defining a dimensionless temperature 〈Ts 〉s -T∞ 〈Tf 〉f -T∞ 〈θs 〉s = , 〈θf 〉f = (4.31) T0 - T∞ T0 - T∞ where T0 is the initial temperature of the system and T∞ is the equilibrium temperature. Similarly, to maintain the same scale for all temperature measurements, the wall temperature gradient is made dimensionless by ΔTw (t) Δθw = (T0 -T∞ ) (4.32). After applying the dimensionless parameters and manipulating the equations to solve for the transient term, the dimensionless solid equation is given by 126 ∂〈θs 〉s ∂2 〈θs 〉s ∂t* = ∂z* - asf Lbed BiLbed <〈θs 〉s -〈θf 〉f = (4.33) where BiLbed is a Biot number with the bed length as the length scale. Boundary conditions are given by ∂〈θs 〉s ∂〈θs 〉s " = " = 0 (4.34). ∂z* z* =0 ∂z* z* =1 The dimensionless fluid equation is given by ∂〈θf 〉f ∂2 〈θf 〉f 1 Lbed ∂〈θf 〉f 2k ∂t* = rα F 2 +ϵ dp 9 Lbed asf Nudp <〈θs 〉s -〈θf 〉f = - Redp Pr ∂z* : - > ϵkw ? raspect rw Δθw G (4.35) ∂z* f with boundary conditions ∂〈θf 〉f 〈θf 〉f ϵD α ? ∂z* - ϵ αbed (ρ* ρin ) ṙ O2 (5.13) bed ∂〈θs 〉s ∂2 〈θs 〉s ΔHRxn L2 ṙ O2 ∂t* = ∂z* - asf LBiL <〈θs 〉s -〈θf 〉f =+ (T (5.14) 0 -T∞ ) keff,s MO2 ∂〈θf 〉f ∂2 〈θf 〉f 1 L ∂〈θf 〉f ∂t* = rα F *2 + ϵ d 9 Lasf Nudp <〈θs 〉s -〈θf 〉f = - Redp Pr ∂z* :G (5.15) ∂z p with boundary conditions ∂ρ * ρ* (z = 0) = 1 , ∂z* J =0 (5.16) z* =1 ∂xO2 xO2 P ? is the reference baseline to calculate the relative change in the oxygen concentration Ar Ref taken from equilibrium when PO2 is approximately 0.2 atm. Integration of Eq. B.12 over the full reduction or oxidation portion of the cycle gives the total reactivity. Error in the total reactivity comes from integrating Eq. B.12 after perturbing the value of v̇ O2 by the uncertainty in the flow rate, done for both the positive and negative perturbations. Error in the drop calorimetry measurements comes from uncertainty in the temperature measurements of the water at the initial and final equilibrium states and from the calibration constant. The calibration constant given by Eq. 2.2 is calculated from the slope of the linear fit of the calibration data collected for alumina. The error in the constant is the standard error in the slope of the fit, calculated in Excel using the LINEST function. Error in the total enthalpy measurements in Figure 2.10 is then calculated using a Kline-McClintock type error analysis derived from Eq. 2.2 as 157 ∂ΔHDrop 2 ∂ΔHDrop 2 μΔH = ±D> ∂C1 ? μ2C + > ∂ΔT ? μ2ΔT (B.13). Drop 1 Error in the total enthalpy measurements then propagates to the energy density calculations (both the total and sensible portion) and comes from the root sum of the squares of the individual error components. Error in the acid solution calorimetry measurements comes from the standard error in the slope for the dissolution experiments and error in the calibration constant, which in turn comes from standard error in the slopes for the dissolution experiments of the calibration materials. All errors are propagated using a Kline-McClintock type error analysis. Error in the calibration constant is given by 2 2 ∂C2 ∂C2 μC = ± KL ΔT M μ2ΔT +L ΔT M μ2ΔT (B.14) 2 ∂( ) ( ) ∂( ) ( ) m MgO m MgO m Mg m Mg where Eq. B.14 is derived from Eq. 2.5 and μ(ΔT) is the standard error in the slope from the m dissolution experiments in Figure B.2. Errors in all enthalpy of dissolution measurements are then calculated by 2 ∂ΔHDiss,i ∂ΔHDiss, i 2 μΔH = ± NO ΔT P μ2ΔT + > ∂C2 ? μ2C (B.15) Diss, i ∂( ) ( ) 2 m i m i where Eq. B.15 is derived from Eq. 2.3 and the subscript i serves as a generic placeholder for the different chemical compounds tested in the calorimeter. Because all standard enthalpy of formation equations (derived using Hess’s Law from the chemical reactions in Table 2.3) are simple addition and subtraction, all errors in the standard enthalpies of formation reported in Figures 2.8 and 2.11 are calculated by 158 ∂ΔH 2 μΔH = ± D∑ni= 1 >∂ΔHf? μ2ΔH (B.16) f i i where the subscript i represents all measured or calculated enthalpy terms in the particular equation for which the error is being calculated. Finally, the error in the energy of chemical reaction calculated by Eq. 2.11 is calculated by 2 μΔH Rxn =± y1 - y2 Dμ2ΔHf, OX + μ2ΔHf, RED (B.17) where ΔHf, OX and ΔHf, RED are the standard enthalpies of formation of the oxidized and reduced magnesium-manganese oxide phases, respectively, from Figure 2.11. B.3. Tabulated Data from Bar Charts Much of the results in Chapter 2 are presented as bar charts. As an alternative form of presentation and in case knowledge of the precise value is required, the data are provided again here in the form of tables. The figure to which each table is connected is listed in the caption. Table B.1. Reactivity. Predicted and measured reactivity for reduction of magnesium- manganese oxide of Mn/Mg molar ratios of 2/3, 1/1, and 2/1 at 1500 °C assuming a reference oxidation equilibrium state of TOX = 1000 °C and PO2,OX = 0.2 atm [29],[36]. Data is presented graphically in Figure 2.1. Units are cm3 of oxygen gas released per gram of fully-oxidized magnesium-manganese oxide solid material. PO2, RED Type 2/1 (cm3O2 g-1) 1/1 (cm3O2g-1) 2/3 (cm3O2g-1) (atm) Theoretical 0.2 30.3 41.5 35.4 Theoretical 0.01 54.1 43.3 35.6 Experimental 0.2 32.3± 2.9 34.3 ± 2.9 34.9 ± 3.4 159 Table B.2. Energy Density. Predicted energy densities for reduction of magnesium- manganese oxide of Mn/Mg molar ratios of 2/3, 1/1, and 2/1 at 1500 °C assuming a reference oxidation equilibrium state of TOX = 1000 °C and PO2,OX = 0.2 atm [29],[36]. Data is presented graphically in Figure 2.2. PO2, RED (atm) 2/1 (kJ kg-1) 1/1 (kJ kg-1) 2/3 (kJ kg-1) 0.2 886 1044 979 0.1 957 1063 979 0.01 1191 1064 979 Table B.3. Acid Solution Calorimeter Validation. Standard enthalpies of formation measured using the acid solution calorimeter, with comparison to the accepted literature values. Data are presented graphically in Figure 2.8. Material Accepted ΔHf (kJ mol-1) [48] Experimental ΔHf (kJ mol-1) |% Error| MnO -384.9 -380.7 ± 2.0 1.1 ZnO -348.3 -340.3 ± 1.3 2.3 Mn2O3 -952.2 -957.3 ± 2.1 0.5 Fe2O3 -825.5 -813.4 ± 0.2 1.5 Table B.4. Drop Calorimetry Results. Total enthalpies of the oxidized and reduced magnesium-manganese oxide phases measured from drop calorimetry, relative to a room- temperature equilibrium state of (T = 25 °C, PO2 = 0.2 atm), for Mn/Mg molar ratio variants of 2/3, 1/1, and 2/1. The reduced phase is quenched from an equilibrium state of (TRED = 1500 °C, PO2,RED = 0.2 atm). The oxidized phase is quenched from an equilibrium state of (TOX = 1000 °C, PO2,OX = 0.2 atm). Data are presented graphically in Figure 2.10. Mn/Mg ΔHOX (kJ kg-1) ΔHRED (kJ kg-1) 2/3 893.6 ± 8.2 1373.3 ± 12.1 1/1 927.1 ± 8.5 1369.9 ± 12.1 2/1 858.4 ± 7.9 1217.2 ± 12.3 160 Table B.5. Acid Solution Calorimetry Results. Standard enthalpies of formation for the oxidized and reduced magnesium-manganese oxide phases measured from acid solution calorimetry for Mn/Mg molar ratio variants of 2/3, 1/1, and 2/1. Data are presented graphically in Figure 2.11. Mn/Mg ΔHf, OX (kJ mol-1MgMnO) ΔHf, RED (kJ mol-1MgMnO) 2/3 1353.4 ± 8.3 1260.2 ± 9.7 1/1 1063.6 ± 7.9 989.5 ± 7.2 2/1 764.2 ± 7.1 696.8 ± 9.7 Table B.6. Mass Change Results. Mass change parameters in the chemical formula for magnesium-manganese oxide measured during the calorimetric process. Data are presented graphically in Figure 2.12. Mn/Mg x y2 (mol O mol-1MgMnO) y1 (mol O mol-1MgMnO) 2/3 1.5 0.5087 0.0173 1/1 1 0.4642 0.0772 2/1 0.5 0.4806 0.1368 Table B.7. Energy Density Results. Energy density measurements for the thermochemical cycle operating between equilibrium states of (TRED = 1500 °C, PO2,RED = 0.2 atm) and (TOX = 1000 °C, PO2,OX = 0.2 atm) for magnesium-manganese oxide with Mn/Mg molar ratios of 2/3, 1/1, and 2/1. Data are presented graphically in Figure 2.13. Mn/Mg ΔHDrop (kJ kg-1) ΔHRxn (kJ kg-1) ΔHS (kJ kg-1) 2/3 479.1 ± 14.6 590.9 ± 62.5 1070 ± 64.2 1/1 442.8 ± 14.8 586.3 ± 55.0 1029 ± 57.0 2/1 358.8 ± 14.6 565.3 ± 54.8 924.1 ± 56.7 161 APPENDIX C: ATTEMPTS AT STEADY STATE INTERSTITIAL HEAT TRANSFER COEFFICIENT MEASUREMENTS This appendix outlines the work done towards measuring the interstitial heat transfer coefficients in packed beds of spheres using a steady-state approach. As mentioned in Chapter 4, the Joule heating method using direct current (DC) was chosen. However, this approach suffered from many complications, and the few results achieved were inconclusive. The work that was done and the results that were acquired are summarized, along with some brief discussion of the observed phenomena that made this method difficult. Finally, some suggestions are made for future avenues of exploration if it is desired to try to make this method more suitable for use. C.1. Equation Derivation & Error Analysis The equations for extracting the interstitial heat transfer coefficients using the DC Joule heating method are derived from the following energy balances. First, an energy balance is taken over the solid particle bed at steady-state. Assuming that all energy in the solid phase is because of Joule heating and is entirely removed by interstitial convection with the fluid, the energy balance is ΔĖ s = 0 =Ṗ DC - Q̇ Int = VDC IDC - hsf asf VCV >〈Ts 〉 s Avg - 〈Tf 〉 f Avg ? (C.1). Here, Ṗ DC is the power applied to the bed, VDC is the DC voltage drop across the bed, IDC is the current, and >〈Ts 〉 s Avg - 〈Tf 〉 f Avg ? is the difference between the average volume-averaged solid and fluid temperatures. Specifically, 〈Ts 〉 s Avg is the solid particle temperature averaged across the entire packed bed, and 〈Tf 〉 f Avg is the fluid temperature averaged across the entire packed bed. Thus, >〈Ts 〉 s Avg - 〈Tf 〉 f Avg ? is simply the difference between these two values. It is inevitable that there will be heat loss to the environment. Thus, not all of the heat generated by the applied power will be captured by the measurement of the average fluid 162 temperature. To quantify the loss, an energy balance over the entire system at steady state is taken. Assuming that all energy within the system is generated by Joule heating of the particles and is either removed by airflow or lost through the side walls, the system-level energy balance is ΔĖ Sys = 0 = Ṗ DC - ΔHf - Q̇ Loss (C.2). Here, ΔHf is the change in the enthalpy of the air across the bed and Q̇ Loss is the overall heat lost from the system to the surrounding. This equation suggests that to correct for heat loss, the enthalpy change of the fluid across the bed must be used as the independent variable in the analysis instead of the DC power input. Assuming that the heat loss occurs only from the solid phase such that the measured fluid temperatures are independent of the heat loss and that the enthalpy change of the fluid across the bed is equal to the net energy transfer from solid to fluid, an energy balance over just the fluid phase at steady state gives ΔĖf = 0 = Q̇ Int - ΔHf = hsf asf VCV >〈Ts 〉 s Avg - 〈Tf 〉 f Avg ? - ṁ f cpf (Tout, Avg - Tin, Avg ) (C.3) where Tout,Avg and Tin,Avg are the average fluid temperature across the face of the outlet and inlet, respectively. Manipulation of Eq. C.3 gives the expression for data extraction. Rearranging to solve for >〈Ts 〉 s Avg - 〈Tf 〉 f Avg ? gives 1 >〈Ts 〉 s Avg - 〈Tf 〉 f Avg ? = >h ? ΔHf (C.4). sf asf VCV Here, hsf is the interstitial heat transfer coefficient, asf is the specific surface of the particle bed calculated from Eq. 3.14, and VCV is the volume of the control volume, which in this case is the entire internal volume of the section of the tube that contains the packed bed. 163 Eq. C.4 takes the form of the equation of a line through the origin, y = mx. If the set of measurements for a constant flow are plotted with ΔHf as the independent variable, a linear fit through the data gives a slope of the fit m containing the HTC, given by 1 mfit = hsf asf VCV (C.5). Defining a Nusselt number to be hsf dp Nudp = kf (C.6) where dp is the particle diameter, rearranging Eq. C.5 to solve for the heat transfer coefficient, and substituting Eq. C.5 into Eq. C.6 gives the expression for extracting the Nusselt number as 1 1 d Nudp = >m ? >a ? > kp ? (C.7). fit sf VCV f Error in the individual datapoints that yield the slope come from the following components: (1) calibration errors in the inlet, outlet, and solid temperature measurements and (2) uncertainty in the flow rate. These data compile into the standard error in the slope, calculated in Excel using the LINEST function. The total error in the Nusselt number can then be calculated as ∂Nud 2 ∂Nudp μNu = ± Nμ2m > ∂m p ? = μm ∂mfit (C.8) dp fit fit fit where μm is the standard error in the slope of the linear fit. The uncertainties associated with fit each measurement are given in Table C.1. 164 Table C.1. Instrumentation Error. Experimental uncertainties and the associated instrument. All temperature probes were calibrated in a water bath to a high-precision Pt-100 thermometer probe, with errors rounded to the nearest higher tenth (i.e., 0.24 round to 0.3) to ensure caution in quantifying the error. Instrument # Device Measurement Error (+/-) 1 T-type thermocouple Solid T, 8mm 0.3 °C 2 T-type thermocouple Solid T, 8mm 0.3°C 3 T-type thermocouple Solid T, 8mm 0.3 °C 4 T-type thermocouple Solid T, 5mm 0.3 °C 5 T-type thermocouple Solid T, 5mm 0.3 °C 6 T-type thermocouple Solid T, 5mm 0.3 °C 7 Pt-100 RTD Inlet T 0.2 °C 8 Pt-100 RTD Inlet T 0.2°C 9 Pt-100 RTD Inlet T 0.2 °C 10 Pt-100 RTD Outlet T 0.2 °C 11 Pt-100 RTD Outlet T 0.2 °C 12 Pt-100 RTD Outlet T 0.2 °C 13 OMEGA 5500A Series Flow rate of air 1.5% of reading C.2. Experimental Setup Steady-state measurements are performed using 304 stainless steel spheres (304SS) of approx. 5 and 8mm in diameter. Stainless steel particles were chosen because preliminary tests showed that the electrical properties of stainless steel at the temperatures of interest (25-75 °C) allowed a reasonably sized bed to be heated with direct current at amperages achievable with a benchtop power supply (10A DC). Particles are packed inside a vertically-oriented 5.1 cm (2”) OD 25% glass-reinforced PTFE tube (ID approx. 4.4 cm), chosen for its high melting temperature (approx. 260 °C) and electrical insulating properties so that both error from heating of the tube wall and the risk of electrical shock from touching the setup are avoided. An inactive 15.2 cm (6”) high bed of 4mm alumina spheres is placed upstream of the test bed to stabilize the flow field and provide support. Two aluminum mesh electrodes are wired to transmit direct current to the bed. The inlet 165 electrode is installed permanently inside the tube and supported by the alumina spheres. A custom electrode part was built for the outlet. The aluminum mesh is bolted inside a thin aluminum tube for support. The outer end of the thin aluminum tube is bolted inside a 0.64 cm (¼”) thick-walled aluminum ring, to which the power supply is wired. All electrical connections are bolted together and coated with a layer of silver conductive epoxy to ensure good connectivity between components. The entire part slides inside the PTFE tube with a tight fit, yet is removeable so that the test bed can be packed through the outlet end of the tube. The electrode part also serves as a packing support mechanism. The mesh inside the electrode is pressed tightly against the top of the bed, holding the bed in place. The top of the electrode is pressed against a PTFE plate with a hole cut through the center. This plate is bolted to a similar aluminum plate that slides over the PTFE tube. Both plates are bolted using moveable brackets fixed to the support structure column built from 1” slotted aluminum framing. The resulting mechanism fixes the bed in place and completes the circuit between the bed and power supply with good electrical contact. A stainless steel compression fitting is attached to the upstream (bottom) end of the setup to flow air through the particle bed. Compressed air from the main building air compressor is regulated at 60psi by a gauge upstream of the setup to remove tank fluctuations. It then passes through a combination oil/particle/water filter setup prior to reaching the flow controller. Flow rates are regulated by an Omega FMA 5500A Series mass flow controller, range 0-50 SLPM. Three Inconel-sheathed, mica insulated Pt-100 RTD probes are placed inside the bed of alumina spheres just upstream of the inlet electrode to measure air temperature across the face of the inlet. 1/8”-NPT stainless steel compression fittings are embedded in the tube wall to seal the probes in place. Three more Pt-100 probes are inserted through the outlet electrode from the top 166 and are fixed to an overhead slotted-framing support frame such that their tips are just downstream of the outlet of the bed. These probes measure the outlet air temperature across the face of the bed. Solid temperatures are measured at heights of 1, 3, and 5 cm along the 6.5cm high bed. To guarantee correct measurement of the solid temperature, holes were drilled in three of the 304SS spheres. Inconel-sheathed, mica-insulated T-type thermocouples are placed inside the holes and affixed using a Loctite epoxy capable of surviving above 100 °C. Prior to installation, all temperature instruments were calibrated to a high-precision Pt-100 thermometer standard within a hot water bath between 22-75 °C. Additional points were taken using ice water and boiling water to expand the calibration range. Temperature data are recorded using a LabJack T7 data acquisition board. Current is measured using a Fluke 45 dual-display bench-top multimeter. The voltage across the bed is measured by crimping thin copper wires to the inlet and outlet mesh electrodes and connecting them across a Fluke-287 multimeter. Power is calculated from the DC current and voltage measurements. Alumina-silica fibrous ceramic insulation is wrapped around the test section to reduce heat loss. A schematic diagram of the experimental setup is shown in Figure C.1. 167 Figure C.1. Experimental Setup. A schematic diagram of the steady state experimental setup is shown. C.3. Experimental Method and Data Extraction The specific method for the steady state measurements is as follows. First, the flow controller is set at the desired flow rate, forcing air to flow through the bed. The power supply is then set at a fixed current. Finally, the DAQ board is set to record the inlet, outlet, and solid temperatures. The system is left undisturbed until the temperature signals reach a steady state. At this point, the average temperature is taken over the previous 100s for each temperature probe. The average solid temperature is calculated from the three solid temperature measurements. The average inlet and outlet temperatures are taken across each face assuming a uniform plug flow. These are then averaged to determine the average fluid temperature. The assumption is that the bed is short enough to assume a linear air temperature profile along the column. Finally, the average phase temperature difference is calculated as the difference between the average solid and average fluid temperatures. These values are plotted as a function of the enthalpy change of the air across the bed. Data are collected at ten different current inputs. A linear fit is taken 168 through the datapoints and the Nusselt number for the given flow rate is extracted from the slope of the fit using Eq. C.6. C.4. Results Results were achieved at only six Reynolds numbers, with flow rates of 25, 35, and 45 SLPM for the 5mm and 8mm sphere beds. The linear fits used to extract the interstitial heat transfer coefficients are shown in Figure C.2. Figure C.2. Datasets. Linear fits of the steady state data as prescribed by Eq. C.4 for the packed beds of 304SS spheres with a diameter of (a) 5mm and (b) 8mm. (a) 12.00 10.00 Re = 92 Re = 170 8.00 (Ts,Avg - Tf,Avg) (°C) Re = 131 6.00 4.00 2.00 0.00 0.00 5.00 10.00 15.00 20.00 25.00 30.00 ∆Hair (W) 169 Figure C.2 (cont’d). (b) 20.00 18.00 Re = 213 16.00 Re = 151 14.00 Re = 89 (Ts,Avg - Tf,Avg) (°C) 12.00 10.00 8.00 6.00 4.00 2.00 0.00 0.00 5.00 10.00 15.00 20.00 ∆Hair (W) Results for the steady state measurements are shown on both Nusselt and Chilton- Colburn factor plots in Figures C.3 and C.4, respectively, with attached error bars. 170 Figure C.3. Results, Nusselt. Nusselt plot of steady state results, with literature studies included for comparison. The limitations of the power supply restricted measurements to Reynolds numbers of approx. 220 and below. The available results are noticeably lower in magnitude than much of the surveyed literature. 60 Wakao/Kagui/Funazkri (1979) 5mm 8mm Gunn (1978) 50 Furnas (1930) Glaser & Thodos (1958) Singhal et al (2016) (DNS Spheres) Achenbach (1995) Gnielinski (1978) 40 Gupta, C haube, Upadhyay (1974) Whitaker (1972) Sen Gupta & Thodos (1962) De Acetis & Thodos (1960) Nudp(-) 30 Malling & Thodos (1966) McConnachie & Thodos (1963) 20 10 0 0 50 100 150 200 250 300 Redp (-) 171 Figure C.4. Results, Chilton-Colburn. Steady state results presented as a Chilton-Colburn factor, with some scanned literature data for comparison. The lower magnitude of the results becomes much more apparent when plotted in this form, as three of the six datapoints do not align at all with the general trend, while the remaining three fall along the lower bounds of the data range. 0.1 !jH (-) Scanned Data Gupta et al (1974) 5mm 8mm 0.01 50 500 5000 Redp(-) It is quite clear from Figures C.3 and C.4 that these results fall on the lower boundary of much of the reported literature. While the results are not entirely out of agreement with the literature, this does raise some concerns as to the validity of the results without much more extensive investigation. Particularly, results at higher Reynolds numbers are necessary before any extensive conclusions can be drawn. C.5. Discussion of Difficulties and Suggestions for Future Work As mentioned, this particular steady state approach had a lot of challenges that raised enough concerns to consider the reported results unreliable. The first such was the lack of a 172 perfectly linear temperature profile along the central axis of the bed. Rather, the three solid temperature measurements indicated that a parabolic profile occurred, with the temperature leveling off along the bed as the outlet was approached. This raises some concerns about the data extraction process because the lumped analysis as-written assumes a linear air temperature profile to rely on just the inlet and outlet temperature measurements. Second, the temperature recordings were inconsistent. Sometimes, there seemed to be a slight radial temperature profile across the face of the outlet, with the coolest temperature at the centerline and a gradient magnitude of approx. 2 °C. At other times, there was a much larger temperature profile where the hottest temperature was on one side of the face, the coolest on the opposite side, and the value at the midpoint of the face in the middle, where the gradient across the face could reach as large as 7-8 °C in the range of datapoints acquired. It is hypothesized that this was because of an uneven flow of current through the particles in the bed, resulting in one portion of the bed experiencing higher heating than the rest. It is possible that this was caused by disruptions in the packing (and thus the circuit) caused by the insertion of the three temperature probes for measuring the solid temperature. In this instance, the average temperature measurements are less reliable, particularly for the solid phase, as the particles with embedded probes were placed along the centerline of the bed. In theory, the use of a linear fit should help balance out any effects like this because it is the slope that is of interest, but it is impossible to say for sure without achieving a uniform current flow and the associated uniform heating. Finally, there were some instances where during the initial heating, the temperature signals of the solid probes would sometimes become very noisy, then suddenly drop to a very stable value. Observations of this effect showed that it was because of the direct current flow, as this happened only when the current was increased and particularly so when the power supply 173 was first turned on. The probes were sheathed and had been fixed to the spheres with an epoxy layer that should have blocked electrical conduction, yet this effect was still observed. Because the signals did eventually stabilize and make sense relative to the other measurements and did not spike when the power supply was turned off under these conditions, it is unclear whether there was a consistent error in the measurement from the current or not, increasing the overall lack of confidence in the system. In general, the concerns listed above can be reduced to the concern that not enough suitably accurate temperature measurements are acquired such that true averages of the solid and fluid temperature fields are measured. Thus, the resulting measurements and the overall method are considered unreliable. There are several possible solutions to these issues that are worth investigating if this method is to be attempted again. First, any setup should be significantly larger than the one used in this work. This would minimize the disruptions to the packing caused by the probes in the bed. Second, probes that are inserted in the bed should be coated with a more substantial layer of an electrically-insulating compound to ensure there is no noise in the measurement signal from the conduction of electricity. Third, more temperature sensors should be placed in the bed to measure the solid phase at multiple locations in both the axial and radial directions in case non- uniform heating still occurs. Finally, probes should be inserted into the pore space as close as possible to the particles in which probes are inserted to try to measure the local temperature difference in the bed. This would remove the reliance on the linear air temperature profile assumption and allow for direct computation of the heat transfer coefficient at many locations, from which an average could be taken. While this certainly does not guarantee success, any future work with this approach should begin by following these suggestions. 174 APPENDIX D: DERIVATION AND SCALING OF THE MODIFIED C-S MODEL This appendix outlines the step-by-step derivation and scaling of the modified C-S model used in Chapter 4 to extract the interstitial heat transfer coefficients from the experimental data and to perform the parametric analysis of the error introduced to the Nusselt number by the transient wall heat transfer effect. D.1. Derivation of the Solid-Phase Energy Equation Define the control volume (CV) as the cross section of the bed only along the length Δz in which the solid phase is contained. The energy balance for the solid phase is given as ∂〈T 〉 s dĖ s =<ρcp = Vs ∂ts = Ė in - Ė out + Ė gen = qz,s - qz+Δz,s -qInt (D.1) s where 𝜌 is the density, cp is the specific heat, Vs is the volume of the control volume occupied by just the solid phase, of the 〈Ts 〉s is the intrinsic volume-averaged solid temperature, qz,s , qz+Δz,s is the axial conductive heat transfer into and out of the control volume due to conduction, respectively, and qInt is the interstitial convection heat transfer with the fluid. Substituting the appropriate area and volume terms into Eq. D.1 gives ∂〈Ts 〉s <ρcp = Vs ∂t = q" z,s Ac - q" z+Δz,s Ac - q'''Int VCV (D.2) s where q" is the axial heat flux, Ac is the cross-sectional area of the control volume, q'''Int is the interstitial convection per unit volume, and VCV is the volume of the control volume. To define both sides of the equation over the same volume, the relationship between the solid volume and the control volume must be established. This is done by defining a generic phase fraction as Vi ϵi = VCV , ∑ni ϵi = 1 (D.3). Because there are only two phases within the control volume, a single phase fraction 𝜖 can be used. Thus, the bulk porosity of the porous medium is defined from Eq. D.3 as 175 Vf Vs Vs ϵ =V =1- VCV →1-ϵ= VCV (D.4). CV Solving Eq. D.4 for Vs and substituting into Eq. D.2 gives ∂〈Ts 〉s <ρcp = (1-ϵ)VCV ∂t = q" z,s Ac - q" z+Δz,s Ac - q'''Int VCV (D.5). s Distributing the volume term to the right-hand side, applying VCV = Ac Δz, and simplifying gives ∂〈Ts 〉s q" z,s Ac q" z+Δz,s Ac q''' VCV (q" z+Δz,s - q"z,s ) <ρcp = (1-ϵ) ∂t = VCV - VCV - Int VCV =- Δz - q'''Int (D.6). s Taking the limit of Eq. D.6 as Δz approaches 0 gives ∂〈Ts 〉s ∂q''z,s <ρcp = (1-ϵ) ∂t =- ∂z - q'''Int (D.7). s Applying Fourier’s Law to the conduction term gives ∂q''z,s ∂ ∂〈Ts 〉s ∂2 〈Ts 〉s - ∂z = - ∂z >-k ∂z ? = keff,s ∂z2 (D.8) where k is defined as the effective thermal conductivity of the bed, keff,s. Applying Newton’s Law of Cooling to the volumetric interstitial heat transfer term gives q'''Int = hV (TH -TC ) = hV <〈Ts 〉s -〈Tf 〉f = (D.9) where hV is defined as the volumetric heat transfer coefficient with units Wm-3 K-1. Substituting Eq. D.8 and Eq. D.9 into Eq. D.7 gives ∂〈Ts 〉s ∂2 〈Ts 〉s <ρcp = (1-ϵ) ∂t = keff ∂z2 - hV <〈Ts 〉s -〈Tf 〉f = (D.10). s To define a heat transfer coefficient per unit surface area instead of per unit volume, the surface area per unit volume, i.e., the specific surface, must be defined. First, the volume taken up by particles is given by Vp Np = VCV (1-ϵ) (D.11) Where Np is the number of particles in a bed and Vp is the volume of a single particle. Rearranging Eq. D.11 gives the number of particles per unit volume as 176 Np 1-ϵ VCV =V (D.12). p Multiplying Eq. D.12 by the surface area of a single particle gives the definition of the specific surface as N A asf = Ap >V p ? = Vp (1-ϵ) (D.13). CV p The interstitial convection heat transfer coefficient per unit surface area is then defined as hV hsf = asf (D.14). Thus, the final form of the solid-phase energy equation is given as ∂〈Ts 〉s ∂2 Ts <ρcp = (1-ϵ) ∂t = keff,s ∂z2 - hV <〈Ts 〉s -〈Tf 〉f = (D.15). s D.2. Derivation of the Fluid-Phase Energy Equation The control volume is again defined as the cross section of the bed along length Δz. The assumption is made that all heat transfer from the wall to the packed bed occurs as convective heat transfer to the fluid phase. The energy balance for the fluid phase is given as ∂〈T 〉 f dĖ f = <ρcp = Vf ∂tf = Ė in - Ė out + Ė gen = qz,f - qz+Δz,f + qInt + qw + ṁ f Hin - ṁ f Hout (D.16) f where qz,f , qz+Δz,f are the axial heat conduction into and out of the control volume, respectively, ṁ f is the mass flow rate of fluid through the control volume, H is the enthalpy of the fluid at the inlet and outlet of the control volume, qw is the convective heat transfer from the wall to the fluid, and qInt is the interstitial convective heat transfer with the solid phase. Substituting the appropriate area and volume terms into Eq. D.16 gives ∂〈Tf 〉f <ρcp = Vf ∂t = q" z,f Ac -q" z+Δz,f Ac + q'''Int VCV + q''w Pw Δz - ṁ f (Hout - Hin ) (D.17) f where q''w is the heat flux across the tube wall and Pw is the perimeter of the tube wall. Solving Eq. D.4 for Vf and substituting in to Eq. D.17 gives 177 ∂〈Tf 〉f <ρcp = ϵVCV ∂t = q" z,f Ac - q" z+Δz,f Ac + q'''Int VCV + q''w Pw Δz - ṁ f (Hout - Hin ) (D.18). f Distributing the volume term, applying VCV = Ac Δz and Pw = 2πR and simplifying gives " " ∂〈Tf 〉f q z,f Ac q z+Δz,f Ac q'''Int VCV q''w (2πRbed )Δz ṁ f <ρcp = ϵ = - + + - (Hout - Hin ) f ∂t VCV VCV VCV VCV VCV (q" z+Δz,s - q''z,s ) 2 q''w - V ḟ (Hout - Hin ) m =- Δz + q'''Int + R (D.19) bed CV where Rbed is the radius of the packed bed contained within the control volume. Defining the mass flux of fluid through the control volume as ṁ m"f = ρuD = A f (D.20), c where uD is the Darcy velocity (i.e., superficial velocity) of the fluid, and the enthalpy change across the control volume as Hout - Hin = ΔH = cp ΔT = cp (〈Tf 〉fout -〈Tf 〉fin ) (D.21), then substituting Eq. D.20 and Eq. D.21 into Eq. D.19 gives ∂〈Tf 〉f (q" z+Δz,s - q''z,s ) " f f 2 q"w mf cp f 4〈Tf 〉out -〈Tf 〉in 5 <ρcp = ϵ ∂t =- Δz + q'''Int + R - Δz (D.22). f bed Simplifying each term and then taking the limit of Eq. D.22 as Δz approaches 0 gives ∂〈Tf 〉f ∂q''z,f 2 ∂〈Tf 〉f <ρcp = ϵ ∂t =- ∂z + q'''Int + R q"w - m"f cpf ∂z (D.23). f bed Applying Fourier’s Law to the conduction term gives ∂q"z,f ∂ ∂〈Tf 〉f ∂2 〈Tf 〉f - ∂z = - ∂z >-k ∂z ? = keff,f ∂z2 (D.24) where keff,f is defined as the fraction of the cross-sectional area occupied by just the fluid phase, 𝜀kf. Applying Newton’s Law of Cooling to the volumetric interstitial heat transfer term gives q'''Int = hV (TH -TC ) = hV <〈Ts 〉s -〈Tf 〉f = (D.25) 178 where hV is again defined as the volumetric interstitial heat transfer coefficient with units Wm- 3 K-1. Applying Fourier’s law to the heat flux across the tube wall gives ∂T k q"w = -kw ∂x w = t w <-ΔTw (t)= (D.26) w w where kw and tw are the thermal conductivity and thickness of the wall, respectively. Here, the assumption is that direct measurement of the heat flux across the tube wall is equivalent to the convective heat transfer from the tube wall to the fluid phase. The spatial derivative of the wall temperature is simplified to an approximation across the wall thickness because (a) the temperature gradient across the tube wall is directly measured and (b) the wall is of a known thickness. Because the wall is very thin, the conduction length is defined as the average wall thickness instead of the natural log of the outer-to-inner radius (ln(rout/rin)), as it is assumed that there is no significant difference between a Cartesian and cylindrical coordinate definition across such a thin length. The temperature gradient across the wall is defined such that a positive temperature gradient is in the direction of positive r (i.e., heat transfer from the bed to the wall), whereas a negative temperature gradient is in the direction of negative r (i.e., heat flowing from the wall to the bed). This is in agreement with the behavior observed by the thermopiles in the experimental setup. Finally, substituting Eq. D.24, Eq. D.25 and Eq. D.26 into Eq. D.23 gives the final form of the fluid-phase energy equation as ∂〈Tf 〉f ∂2 〈Tf 〉f 2 kw ∂〈Tf 〉f <ρcp = ϵ ∂t = ϵkf ∂z2 + hV <〈Ts 〉s -〈Tf 〉f =+ R <-ΔTw (t)= - m"f cp (D.27). f bed tw f ∂z D.3. Boundary and Initial Conditions The boundary conditions for the solid equation assume zero heat flux across both the inlet and the outlet boundaries. The boundary conditions for the fluid equation prescribe a known temperature at the inlet and assume zero heat flux across the outlet boundary. The initial 179 condition is assumed to be a uniform temperature and that both phases begin in thermal equilibrium. These are given mathematically as ∂〈Ts 〉s ∂〈Ts 〉s -keff,s ∂z @ = keff,s ∂z @ =0 (D.28) z=0 z =L ∂〈Tf 〉f 〈Tf 〉f (z=0) = Tinlet , -kf J =0 (D.29) ∂z z=L 〈Ts 〉s (z, t = 0) = 〈Tf 〉f (z, t = 0) = T0 (D.30) where Lbed is the length of the packed bed. D.4. Scaling the Modified C-S Model To make the modified C-S model dimensionless, a set of scaling parameters must be defined. For the time scale, a parameter similar to a Fourier number based on the properties of the bed as keff,s αbed L2 t* = t= t→ t = αbed t* (D.31) (1-ϵ)ρs cp L2bed L2bed bed s where 𝛼bed is the thermal diffusivity of the packed bed, defined as keff,s αbed = (1-ϵ)ρcp (D.32). The dimensionless length scale is defined as z z* = L → z = Lbed z* (D.33). bed The dimensionless volume-averaged temperature for phase i is defined as 〈Ti 〉i - T∞ 〈θi 〉i = (D.34) T0 - T∞ where T0 and T∞ are the initial and final temperatures, respectively. Substituting Eq. D.31, Eq. D.33, and Eq. D.34 into Eq. D.15 gives 180 ∂(〈θs 〉s (T0 -T∞ )+T∞ ) ∂2 (〈θs 〉s (T0 -T∞ )+T∞ ) (1-ϵ)ρs cps =k eff,s (1-ϵ)ρs cp L2bed ∂(z* Lbed )2 ∂L s tM* keff,s −hV <(〈θs 〉s (T0 -T∞ )+T∞ )-(〈θf 〉f (T0 -T∞ )+T∞ )= (D.35) where 〈θs 〉s is the dimensionless intrinsic volume-averaged solid temperature. Removing the constants from the derivatives, deleting the derivatives of constants, and simplifying gives keff,s ∂〈θs 〉s keff,s ∂2 〈θs 〉s = - hV <〈θs 〉s -〈θf 〉f = (D.36). L2bed ∂t* L2bed ∂z* Rearranging Eq. D.36 to isolate the transient term gives ∂〈θs 〉s ∂2 〈θs 〉s hV L2bed ∂t* = ∂z* - k <〈θs 〉s -〈θf 〉f = (D.37). eff,s Defining a bed-based Biot number as hsf Lbed BiLbed = keff,s (D.38) and substituting Eq. D.38 into Eq. D.37 gives the final dimensionless solid phase energy equation as ∂〈θs 〉s ∂2 〈θs 〉s ∂t* = ∂z* - asf Lbed BiLbed <〈θs 〉s -〈θf 〉f = (D.39) with dimensionless boundary conditions ∂〈θs 〉s ∂〈θs 〉s ∂z* @* = ∂z* @* =0 (D.40). z =0 z =1 The scaling parameters for the fluid phase energy equation must match solid equation for both to be on the same scale. Therefore, substituting Eq. D.31, Eq. D.33, and Eq. D.24 into Eq. D.27 gives 181 ∂(〈θf 〉f (T0 -T∞ )+T∞ ) ∂2 (〈θf 〉f (T0 -T∞ )+T∞ ) ϵρf cpf (1-ϵ)ρs cp L2bed = ϵkf 2 + hV <(〈θs 〉s (T0 -T∞ )+T∞ -(〈θf 〉f (T0 -T∞ )+T∞ )= s ∂4z* Lbed 5 ∂E t* F keff,s ∂(〈θf 〉f (T0 -T∞ )+T∞ ) 2 kw - m"f cpf ∂4z* Lbed 5 +R <-ΔTw (t)= bed tw (D.41) where 〈θf 〉f is the dimensionless intrinsic volume-averaged fluid temperature. Removing the constants from the derivatives, deleting the derivatives of constants, and simplifying gives 2 k w ϵρf cp keff,s ∂〈θf 〉f ϵkf ∂2 〈θf 〉f m"f cp ∂〈θf 〉f Rbed tw 4-ΔTw (t)5 f = + hV <〈θs 〉s 〉f -〈θf =- L f + (D.42). (1-ϵ)ρs cp L2bed ∂t* L2bed ∂z* 2 bed ∂z * (T0 -T∞ ) s L2bed Multiplying all terms in Eq. D.42 by ϵkf gives k w ρf cp keff,s ∂〈θf 〉f ∂2 〈θf 〉f hV L2bed m"f cp Lbed ∂〈θf 〉f 2L2bed tw 4-ΔTw (t)5 (1-ϵ)ρs cp f kf ∂t* = 2 + ϵkf <〈θs 〉s -〈θf 〉f =- f ϵkf ∂z* + ϵR k (T -T ) (D.43). s ∂z* bed f 0 ∞ Rearranging Eq. D.43 to isolate the transient term gives ∂〈θf 〉f (1-ϵ)ρs cp kf ∂2 〈θf 〉f hV L2bed m"f cp Lbed ∂〈θf 〉f 2L2bed kw 4-ΔTw (t)5 ∂t* = ρf cp s keff,s F *2 + ϵkf <〈θs 〉s -〈θf 〉f =- f ϵkf ∂z* + ϵR (T0 -T∞ ) G (D.44). f ∂z bed kf tw To simplify Eq. D.44 into a collection of dimensionless parameters, the ratio of fluid-to-bed- thermal diffusivities is defined as (1-ϵ)ρs cp kf αf s rα = ρf cp keff,s =α (D.45), f bed a Nusselt number based on bed length and a Nusselt number based on particle diameter as hsf Lbed NuLbed = kf (D.46) hsf dp Nudp = kf (D.47), and a Reynolds number based on bed length and a Reynolds number based on particle diameter as 182 m"f Lbed cpf μf ReLbed Pr = μf kf (D.48) m"f dp cpf μf Redp Pr = μf kf (D.49) where Pr is the Prandtl number and appears by incorporating the dynamic viscosity 𝜇f and assists in making the equation dimensionless. Relating Eq. D.46 to Eq. D.47 and Eq. D.48 to Eq. D.49 through a change of length scale, then substituting the resulting terms and Eq. D.45 into Eq. D.44 gives ∂〈θf 〉f ∂2 〈θf 〉f asf Lbed Lbed 1 Lbed ∂〈θf 〉f 2L2bed kw 4-ΔTw (t)5 ∂t* = rα F 2 + Nudp <〈θs 〉s -〈θf 〉f =- ϵ Redp Pr ∂z* + ϵR G (D.50). ∂z* ϵ dp dp bed kf tw (T0 -T∞ ) Simplifying Eq. D.50 by regrouping terms and defining two dimensionless geometry ratios as Lbed raspect = Rbed (D.51) Lbed rw = tw (D.52) and defining a dimensionless wall temperature gradient as ΔTw (t) Δθw = (T0 -T∞ ) (D.53) gives the final dimensionless form of the fluid phase energy equation as ∂〈θf 〉f ∂2 〈θf 〉f 1 Lbed ∂〈θf 〉f 2k ∂t* = rα F *2 +ϵ dp 9 Lbed asf Nudp <〈θs 〉s -〈θf 〉f = - Redp Pr ∂z* : - > ϵkw? raspect rw Δθw G (D.54) ∂z f with dimensionless boundary conditions ∂〈θf 〉f 〈θf 〉f ϵkw? raspect rw Δθw G (D.54) ∂z* f ∂〈θf 〉f 〈θf 〉f