NUMERICAL STUDY OF THE EFFECTS OF FINS AND THERMAL FLUID VELOCITIES ON THE STORAGE CHARACTERISTICS OF A CYLINDRICAL LATENT HEAT ENERGY STORAGE SYSTEM by Wilson Ogoh Submitted in partial fulfillment of the requirements for the degree of Master of Applied Science at Dalhousie University Halifax, Nova Scotia July 2010 © Copyright by Wilson Ogoh, 2010 ii DALHOUSIE UNIVERSITY DEPARTMENT OF MECHANICAL ENGINEERING The undersigned hereby certify that they have read and recommend to the Faculty of Graduate Studies for acceptance a thesis entitled “Numerical Study of the Effects of Fins and Thermal Fluid Velocities on the Storage Characteristics of a Cylindrical Latent Heat Energy Storage System” by Wilson Ogoh in partial fulfillment of the requirements for the degree of Master of Applied Science. Dated: July 27, 2010 Supervisor: _________________________________ Readers: _________________________________ _________________________________ iii DALHOUSIE UNIVERSITY DATE: July 27, 2010 AUTHOR: Wilson Ogoh TITLE: Numerical Study of the Effects of Fins and Thermal Fluid Velocities on the Storage Characteristics of a Cylindrical Latent Heat Energy Storage System DEPARTMENT OR SCHOOL: Department of Mechanical Engineering DEGREE: MASc CONVOCATION: October YEAR: 2010 Permission is herewith granted to Dalhousie University to circulate and to have copied for non-commercial purposes, at its discretion, the above title upon the request of individuals or institutions. _______________________________ Signature of Author The author reserves other publication rights, and neither the thesis nor extensive extracts from it may be printed or otherwise reproduced without the author’s written permission. The author attests that permission has been obtained for the use of any copyrighted material appearing in the thesis (other than the brief excerpts requiring only proper acknowledgement in scholarly writing), and that all such use is clearly acknowledged. iv TABLE OF CONTENTS LIST OF TABLES ............................................................................................................. vi LIST OF FIGURES ......................................................................................................... viii ABSTRACT ..................................................................................................................... xiv LIST OF SYMBOLS USED ............................................................................................. xv ACKNOWLEDGEMENTS ............................................................................................ xvii CHAPTER 1: INTRODUCTION ....................................................................................... 1 1.1 LHESS Application .................................................................................................. 2 1.2 Literature Review ..................................................................................................... 3 1.2.1 Thermal Energy Storage .................................................................................... 3 1.2.2 Phase Change Materials..................................................................................... 5 1.2.3 Melting Phase Change Heat Transfer ................................................................ 9 1.3 Problem Statement ................................................................................................. 15 1.4 Scope of the Thesis ................................................................................................ 17 CHAPTER 2: MODEL GEOMETRY AND MATHEMATICAL FORMULATION ..... 18 2.1 LHESS Geometry and Materials ............................................................................ 18 2.2.1 Fluid Flow Process .......................................................................................... 20 2.2.2 Heat Transfer Process ...................................................................................... 21 2.2.3 Phase change heat transfer ............................................................................... 24 2.3 Numerical Method.................................................................................................. 25 2.3.1 Numerical Modeling ........................................................................................ 25 CHAPTER 3: VALIDATION .......................................................................................... 30 3.1 Fluid Flow Validation ............................................................................................ 30 3.1.1 Analytical Solution .......................................................................................... 30 3.1.2 Numerical Analysis and Validation ................................................................. 32 3.2 Heat Transfer Validation ....................................................................................... 38 3.2.1 Analytical Solution .......................................................................................... 38 3.2.3 Validation ........................................................................................................ 50 3.3.1 Analytical Solution .......................................................................................... 55 v 3.3.2 Numerical Analysis ......................................................................................... 57 3.3.3 Validation ........................................................................................................ 62 CHAPTER 4: LHESS NUMERICAL ANALYSIS: RESULTS AND DISCUSSION .... 66 4.1 Calculation of Energy Stored in the LHESS .......................................................... 66 4.2 Effects of Thermal Fluid Velocity ......................................................................... 67 4.3 Effects of addition of Fins ...................................................................................... 81 4.4 Energy Stored versus Maximum Storage Capacity and Energy Storage System Efficiency ...................................................................................................................... 89 CHAPTER 5: CONCLUSION ......................................................................................... 92 5.1 Conclusion .............................................................................................................. 92 5.2 Recommendations .................................................................................................. 94 REFERENCES ................................................................................................................. 96 APPENDIX A: Temperature Plots as a Function of the Thermal Fluid Velocity .......... 102 APPENDIX B: Temperature Plots as a Function of Number of Fins (HTF Velocity of 0.05m/s) .......................................................................................................................... 109 APPENDIX C: Temperature Plots as a Function of Number of Fins (HTF Velocity of 0.5m/s) ............................................................................................................................ 112 vi LIST OF TABLES Table 1.1: Important characteristics of energy storage materials (Zalba B. et al., 2003) ... 3 Table 2.1: Materials and Dimensions of the LHESS Geometry ....................................... 19 Table 2.2: Thermophysical Properties of Materials .......................................................... 20 Table 2.3: The PCM Fluid Properties ............................................................................... 23 Table 3.1: Initial and boundary conditions for the numerical fluid flow model. .............. 33 Table 3.2: Element sizes, number of elements, DOF and simulation time for the numerical fluid flow model. ...................................................................................... 33 Table 3.3: Initial and boundary conditions for the heat transfer study. ............................ 40 Table 3.4: Element sizes, number of elements, DOF and simulation time for the numerical heat transfer model. .................................................................................. 42 Table 3.5: Fluid reference temperatures cross section B-B, C-C and D-D. ...................... 48 Table 3.6: Pipe wall temperatures at point x, y and z. ...................................................... 48 Table 3.7: Heat fluxes at points x, y and z. ...................................................................... 48 Table 3.8: Numerical convection coefficients (h) ............................................................. 49 Table 3.9: Calculated parameters using turbulent convection in pipe approach. Simulated Time: 5 hours. ............................................................................................................ 50 Table 3.10: Calculated parameters using turbulent convection in pipe approach. Simulated Time: 10 hours.......................................................................................... 50 Table 3.11: Calculated parameters using turbulent convection in pipe approach. Simulated Time: 15 hours.......................................................................................... 50 Table 3.12: Calculated parameters using turbulent convection in pipe approach. Simulated Time: 20 hours.......................................................................................... 51 Table 3.13: Analytical convection coefficients using turbulent convection in pipe approach. .................................................................................................................... 51 Table 3.14: Calculated analytical parameters using turbulent flow over plate with uniform surface heat flux approach. Simulated Time: 5 hours. .............................................. 52 Table 3.15: Calculated analytical parameters using turbulent flow over plate with uniform surface heat flux approach. Simulated Time: 10 hours. ............................................ 53 vii Table 3.16: Calculated analytical parameters using turbulent flow over plate with uniform surface heat flux approach. Simulated Time: 15 hours. ............................................ 53 Table 3.17: Calculated analytical parameters using turbulent flow over plate with uniform surface heat flux approach. Simulated Time: 20 hours. ............................................ 53 Table 3.18: Analytical convection coefficients using turbulent flow over plate with uniform surface heat flux approach. ................................................................................. 53 Table 4.1: Percentage total stored energy increase in the LHESS as thermal fluid velocity is increase from 0.05 m/s to 1 m/s for various numbers of fins. ............................... 70 Table 4.2: Percentage sensible energy increase in the LHESS as thermal fluid velocity is increase from 0.05 m/s to 1 m/s for various numbers of fins. ................................... 76 Table 4.3: Percentage latent energy increase in the LHESS as thermal fluid velocity is increase from 0.05 m/s to 1 m/s for various numbers of fins. ................................... 78 Table 4.4: Maximum stored sensible, latent and total energy in the LHESS as a function of the number of fins. ................................................................................................ 87 viii LIST OF FIGURES Figure 1.1: Schematic representation of LHESS used for hot water storage...................... 2 Figure 1.2: Classification of phase change materials.......................................................... 6 Figure 2.1: (a) Geometry of cylindrical LHESS; (b) Some examples of fin configurations in the cylindrical LHESS (0 to 3 fin arrangement are shown) .................................. 19 Figure 2.2: Modified heat capacity of paraffin wax as a function of temperature. ........... 29 Figure 3.1: Schematic representation of Flow between two plates .................................. 30 Figure 3.2: Numerical model geometry and finite element mesh used to solve for a flow of water between two parallel plates. ........................................................................ 32 Figure 3.3: Position of the presented and analyzed velocity profiles. .............................. 34 Figure 3.4: Simulated fluid velocity profile, for various element sizes, taken on section h- h. Simulated time: 28 hours. ...................................................................................... 34 Figure 3.5: Maximum fluid velocity as a function of the element size used for the numerical fluid flow model at length 49 m. Simulated Time: 28 hours. ................... 35 Figure 3.6: Simulated fluid velocity profiles at different section. Simulated time: 28 hours. ......................................................................................................................... 36 Figure 3.7: Comparison of the analytical and numerical velocity profiles. ...................... 37 Figure 3.8: Numerical model geometry and finite element mesh used to solve convection heat transfer in a pipe. ................................................................................................ 41 Figure 3.9: Simulated temperature plot for the numerical model. Simulated time: 20 hours. ......................................................................................................................... 42 Figure 3.10: Temperature profile, for various element sizes, at section A-A. Simulated Time: 20 hours. .......................................................................................................... 43 ix Figure 3.11: Minimum temperature as a function of the element size, at section A-A. Simulated time: 20 hours. .......................................................................................... 44 Figure 3.12: Simulated thermal fluid velocity profiles at sections B-B, C-C and D-D. Simulated time: 20 hours ........................................................................................... 45 Figure 3.13: Simulated thermal fluid temperature profiles at sections B-B, C-C and D-D. Simulated Time: 20 hours.......................................................................................... 46 Figure 3.14: Schematic thermal fluid temperature profile along sections B-B, C-C and D- D. ............................................................................................................................... 46 Figure 3.15: Numerical convection coefficients as a function of time at sections B-B, C-C and D-D. .................................................................................................................... 49 Figure 3.16: Analytical (turbulent convection in pipe approach) and numerical convection coefficients at different sections as a function of time. ............................................. 52 Figure 3.17: Analytical (turbulent flow over plate with uniform surface heat flux approach) and numerical convection coefficients at different sections as a function of time. ........................................................................................................................... 54 Figure 3.18: Schematic of the 1D analytical phase change problem studied. .................. 55 Figure 3.19: Schematic of the 2D numerical model used to validate phase change heat transfer. ...................................................................................................................... 58 Figure 3.20: Simulated temperature plot for the numerical model. Simulated Time: 21 hours .......................................................................................................................... 59 Figure 3.21: Temperature profile for various element sizes, along section P-P. Simulated Time: 3 hours. ............................................................................................................ 60 Figure 3.22: Temperature distribution as a function of the distance from the wall, along section P-P, for various simulated times. .................................................................. 61 Figure 3.23: Temperature as a function of time, at point B, for a PCM with a melting temperature range of 313K-316K. Simulated Time: 21 hours. ................................. 62 Figure 3.24: Numerical and analytical temperature profile in the PCM along section P-P for various time. ......................................................................................................... 63 Figure 3.25: Numerical and analytical temperature profile along section P-P for various numerical melting temperature ranges. Simulated Time: 21 hours. .......................... 64 x Figure 3.26: Numerical and analytical melting front position along section P-P as a function of time ......................................................................................................... 65 Figure 4.1: Simulated temperature distribution plots for the 6 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. ...................................... 68 Figure 4.2: Simulated temperature distribution plots for the 15 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. ................................ 69 Figure 4.3: Total energy stored in the 0 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s ............................................ 71 Figure 4.4: Total energy stored in the 1 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. ........................................... 71 Figure 4.5: Total energy stored in the 2 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. ........................................... 72 Figure 4.6: Total energy stored in the 3 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. ........................................... 72 Figure 4.7: Total energy stored in the 4 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. ........................................... 73 Figure 4.8: Total energy stored in the 5 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. ........................................... 73 Figure 4.9: Total energy stored in the 6 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. ........................................... 74 Figure 4.10: Total energy stored in the 9 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s....................................... 74 Figure 4.11: Total energy stored in the 12 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s....................................... 75 Figure 4.12: Total energy stored in the 15 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s....................................... 75 Figure 4.13: Total energy stored in the 18 fin LHESS arrangement as a function of the time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1m/s. .............................. 76 Figure 4.14: Sensible heat energy stored after 12 hours in 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin LHESS arrangements for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. ............................................................................................................................ 77 xi Figure 4.15: Latent heat energy stored after 12 hours in 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin LHESS arrangements for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. ................................................................................................................................... 78 Figure 4.16: Total energy stored after 12 hours in 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin LHESS arrangements for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1m/s. ... 79 Figure 4.17: Melting fraction as a function of the number of fins for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. ................................................................. 80 Figure 4.18: Simulated temperature distribution plots for the 1, 5, 10, 14 and 18 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.05 m/s. ................................................................................................................................... 82 Figure 4.19: Total energy stored in some of the 0 to 27 fin LHESS arrangements as a function of time. Thermal fluid velocity of 0.05 m/s. ............................................... 83 Figure 4.20: Sensible and latent heat energy stored in the 0 to 27 fin LHESS arrangement after 12 hours of charging. Thermal fluid velocity of 0.05 m/s................................. 84 Figure 4.21: Simulated temperature distribution plots for the 1, 5, 10, 14 and 18 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.5 m/s. ................................................................................................................................... 85 Figure 4.22: Total energy stored in some of the 0 to 27 fin LHESS arrangements as a function of time. Thermal fluid velocity of 0.5 m/s. ................................................. 86 Figure 4.23: Sensible and Latent heat energy stored in the 0 to 27 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.5 m/s. ........... 88 Figure 4.24: Total energy stored in some of the 0 to 27 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.5 m/s. ............................................... 88 Figure 4.25: Comparison of total energy stored with maximum storage capacity for the 0 to 27 fin LHESS arrangement after 12 hours of charging and with thermal fluid velocities of 0.05 and 0.5 m/s. ................................................................................... 90 Figure 4.26: LHESS storage efficiencies for 0 to 27 fins LHESS arrangements after 12 hours of charging with thermal fluid velocities of 0.05 and 0.5 m/s. ........................ 91 Figure A.1: Simulated temperature distribution plots for the 0 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .................................... 102 xii Figure A.2: Simulated temperature distribution plots for the 1 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .................................... 103 Figure A.3: Simulated temperature distribution plots for the 2 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .................................... 103 Figure A.4: Simulated temperature distribution plots for the 3 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .................................... 104 Figure A.5: Simulated temperature distribution plots for the 4 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .................................... 104 Figure A.6: Simulated temperature distribution plots for the 5 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .................................... 105 Figure A.7: Simulated temperature distribution plots for the 6 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .................................... 105 Figure A.8: Simulated temperature distribution plots for the 9 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .................................... 106 Figure A.9: Simulated temperature distribution plots for the 12 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .............................. 106 Figure A.10: Simulated temperature distribution plots for the 15 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .............................. 107 Figure A.11: Simulated temperature distribution plots for the 16 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .............................. 107 Figure A.12: Simulated temperature distribution plots for the 17 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .............................. 108 Figure A.13: Simulated temperature distribution plots for the 18 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. .............................. 108 Figure B.1: Simulated temperature distribution plots for 0, 1, 2, 3 and 4 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.05 m/s. ....... 109 Figure B.2: Simulated temperature distribution plots for 5, 6, 7, 8 and 9 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.05 m/s. ....... 110 Figure B.3: Simulated temperature distribution plots for 10, 11, 12, 13 and 14 fin LHESS arrangements after 12 hours of charging. Thermal fluid vel. of 0.05 m/s. .............. 110 xiii Figure B.4: Simulated temperature distribution plots for 15, 16, 17 and 18 fin LHESS arrangements after 12 hours of charging. Thermal fluid vel. of 0.05 m/s. .............. 111 Figure C.1: Simulated temperature distribution plots for 0, 1, 2, 3 and 4 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.5 m/s. ......... 112 Figure C.2: Simulated temperature distribution plots for 5, 6, 7, 8 and 9 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.5 m/s. ......... 113 Figure C.3: Simulated temperature distribution plots for 10, 11, 12, 13 and 14 fin LHESS arrangements after 12 hours of charging. Thermal fluid vel. of 0.5 m/s. ................ 113 Figure C.4: Simulated temperature distribution plots for 15, 16, 17 and 18 fin LHESS arrangements after 12 hours of charging. Thermal fluid vel. of 0.5 m/s. ................ 114 xiv ABSTRACT This thesis work presents a numerical study of the effects of fins and thermal fluid velocities on the storage characteristics of a cylindrical latent heat energy storage system (LHESS). The work consists of two main components: 1. The development of a numerical method to study and solve the phase change heat transfer problems encountered in a LHESS during charging of the system, which results in melting of the phase change material (PCM). The numerical model is based on the finite element method. The commercial software COMSOL Multiphysics was used to implement it. The effective heat capacity method was applied in order to account for the large amount of latent energy stored during melting of a PCM, and the moving interface between the solid and liquid phases. The fluid flow, heat transfer and phase change processes were all validated using known analytical solutions or correlations. 2. Due to the low thermal conductivity of PCMs, the heat transfer characteristics of an enhanced LHESS was studied numerically. The effects of fins and the thermal fluid velocity on the melting rate of the PCM in the LHESS were analyzed. Results obtained for configurations having between 0 and 27 fins show that the heat transfer rate increases with addition of fins and thermal fluid velocity. The effect of the HTF velocity was observed to be small with few fin configurations since the thermal resistance offered by the LHESS system, mostly PCM, is vastly more important under these conditions; while its effect becomes more pronounced with addition of fins, since the overall thermal resistance decreases greatly with the addition of fins. The total energy stored after 12 hours for 0 and 27 fins configurations range between 3.6 MJ and 39.7 MJ for a thermal fluid velocity of 0.05 m/s and between 3.7 MJ and 57 MJ for a thermal fluid velocity of 0.5 m/s. The highest system efficiencies for the 0.05 m/s and 0.5 m/s, obtained with 27 fins configuration are 68.9% and 97.9% respectively. xv LIST OF SYMBOLS USED Dimensional variables Cp Heat capacity (J kg -1 K -1 ) D Diameter (m) E Total energy stored (J) g Acceleration of gravity (m s -2 ) h Thickness (m) h Heat transfer coefficient (W m -2 K -1 ) k Thermal conductivity (W m -1 K -1 ) L Latent heat of fusion (J kg -1 ) m Mass (kg) P Pressure (Pa) Q Energy stored (J) r Radius (m) r Cylindrical coordinate (m) T Temperature (K) t Time (s) u Thermal fluid velocity (m s -1 ) u Velocity of the melting front (m s -1 ) u x velocity component (m s -1 ) u HTF velocity (m s -1 ) V Volume (m 3 ) X Solid- liquid interface position (m) x Coordinate (m) y Coordinate (m) z Cylindrical coordinate (m) Greek symbols  Thermal diffusivity of melt (m2 s-1) he Endothermic heat of reaction (J kg -1 ) T Temperature difference (K)  Viscous dissipation (m2 s-2)  Dynamic viscosity (N s m-2)  Kinematic viscosity (m2 s-1)  Density (kg m-3) Non-dimensional variables f Friction factor n Normal direction xvi Subscripts ap Average specific heat eff Effective f Final i Initial l Latent l Liquid in Inlet m Melting max Maximum min Minimum pcm Phase change material s Sensible s Solid T Total w Wall x in the x direction z in the z direction 1 Inner radius 2 Outer radius Definitions of non-dimensional variables Nu Nusselt number ( ) Pr Prandtl number ( ) Re Reynold number ( ) Ste Stefan number ( ) Ra Rayleigh number xvii ACKNOWLEDGEMENTS I thank God for His great goodness, mercy and favours upon me. I would also like to express my sincere gratitude to my supervisor, Professor Dominic Groulx, for his continuous guidance, valuable advice and supervision during the development of this thesis research. He is always ready to help at all times. He is a dedicated supervisor. My special thanks also go to my examining committee members, Professor Peter Allen and Professor Noureddine Ben-Abdallah for their suggestions and comments. I also appreciated the Department of Mechanical Engineering for their educational support throughout my studies. My special thanks are extended to all my friends and colleagues in the Department. I want to use this opportunity to thank my family members for their financial support throughout this period of my studies. I also would like to acknowledge the financial support from the Dalhousie Faculty of Graduate Studies and the Bruce and Dorothy Rosetti Engineering Research Scholarship (2008). Finally, I am greatly indebted to my wife, Elizabeth and my children, Ebenezer, Deborah and Ezra for their understanding, patience and encouragement. I love you and God bless you all. 1 CHAPTER 1: INTRODUCTION During the past twenty years, solid-liquid phase change heat transfer problems associated with latent heat energy storage system (LHESS) have gained much attention because of their extensive usage in industrial refrigeration (Go et al., 2004), crystal growth (Zalba et al., 2003), solar energy units (Hasnain, 1998), welding and casting (Atul et al. 2009). This has also been a popular research area because of the industrial and domestic applications, such as energy recovery from air-conditioning systems (Go et al., 2004), under floor electric heating (Liu et al., 2005), greenhouse heating and other applications (Demirel et al., 1993). The heat capacity of a substance increases by more than one hundred times during solid- liquid phase change transition (Esam et al., 2004); therefore a substantial amount of energy can be stored in the substance resulting from this highly increased thermal storage capacity. The substance used is called a phase change material (PCM). They are widely used in storage and control systems. During low energy demand periods, energy is stored and released during high energy demand periods (Rosen, 2003). Several experimental and theoretical investigations were devoted to modeling the thermal performance of storage units using PCM, solid-liquid phase change heat transfer, investigating new geometries, new concepts for energy storage and transfer and upgrading the available technology (Agyenim et al., 2008). Latent heat thermal storage has been prove to be an effective method for utilization of clean energy sources, because of its high storage density and small temperature variation from storage to extraction (Atul et al. 2009). Energy loss during heat storage and recovery is smaller when compared with that of sensible heat storage; this is due to the finite differences between fluid temperature and PCM (Rosen, 2003). 2 1.1 LHESS Application The application envisioned for the LHESS modeled and studied in this work is thermal storage in a solar hot water system presented in Fig. 1.1. Most solar energy systems require thermal energy storage to eliminate the mismatch between energy supply and demand (Agyenim et al., 2009). LHESS promises high performance and reliability with the advantages of high storage density and nearly constant temperature energy delivery. During the day, when solar energy is available, cold water supplied from the water tank is pumped to the solar collector (valve leading to the LHESS closed). The heated water from the solar collector charges (melts the PCM) the LHESS by flowing through it. The discharging (solidification) process, although not the focus of this research work, takes place at night when solar energy is not available, with the valve leading to the solar collector closed, the cold water pumped from the water tank flows through the LHESS removing the stored thermal energy from the LHESS. Figure 1.1: Schematic representation of LHESS used for hot water storage Water Tank LHES S Supply To taps 3 1.2 Literature Review 1.2.1 Thermal Energy Storage Surveys conducted by several researchers show that thermal energy storage (TES) is one of the most important energy technologies and it has been brought to the forefront of energy research because of its economic benefits stemming from its design and operation in energy conversion systems (Rosen et al., 2003; Anica, 2005; Nayak et al., 2006; Atul et al. 2009; Hasnain, 1998; Zhang et al., 1996). TES systems can collect energy in order to shift its delivery to a later time, or to smooth out the plant output during intermittently cloudy weather conditions. Times of mismatch between energy supply by the sun and energy demand can be reduced, which stand as a key component of any successful thermal system (Agyenim et al., 2008). Thermal energy storage has the potential to produce significant benefits and savings particularly for low-temperature heating and cooling applications. These benefits should see thermal energy storage gain wider acceptance (Dincer, 2002). Good thermal energy storage should allow for minimum thermal energy losses leading to energy savings, while permitting the highest possible extraction efficiency of the stored thermal energy (Hasnain, 1998). Table 1.1 illustrates the important characteristics of energy storage materials based on thermal, physical, chemical and economical properties. These characteristics form the main criteria that govern the selection of phase change heat storage materials. Table 1.1: Important characteristics of energy storage materials (Zalba et al., 2003) Thermal Properties Physical properties Chemical properties Economic properties Phase change temperature fitted to application. High change of enthalpy near temperature of use Low density variation High density. Small or no under cooling. Stability No phase separation Compatibility with container materials. Non-toxic Non-flammable Non- polluting. Inexpensive Abundant 4 Three types of TES can usually be employed: sensible, thermochemical and latent heat storage. Sensible Heat Storage In sensible heat storage (SHS), thermal energy is stored by raising the temperature of a substance (solid or liquid). A sensible thermal energy storage system consists of a storage medium (e.g., water, oil, bricks, sand, soil, or rock beds), a container and input/output devices. The amount of heat stored depends on the specific heat of the medium, the temperature change and the amount of storage material. Therefore, the amount of energy stored Q in a mass m of material is given by: Where and are the final and initial temperatures (K), is the temperature dependent specific heat of the material (J/kg∙K) and is the average specific heat between and (J/kg∙K). Sensible thermal energy storage is mainly used for low- grade heat such as waste heat from power generation plants and industrial thermal processes (Jegadheeswaran et al., 2009), for short and long term storage purposes. Thermochemical heat storage Thermochemical systems rely on the ability of certain materials to absorb and release energy by breaking and reforming bonds at the molecular level in a completely reversible chemical reaction. In this case, the energy stored depends on the amount of storage material, the endothermic heat of reaction, and the extent of conversion, illustrated by the following equation: Where is the fraction reacted, is the mass of heat storage medium (kg) and is the endothermic heat of reaction (J/kg). 5 Although thermochemical systems are not currently viable, a variety of reactions are being explored. At present, lower-temperature reactions (<573K) do not appear promising (Dutt et al., 2004). Latent Heat Storage (LHS) The internal energy associated with the phase of a system during a phase change process is called the latent heat or latent energy. Latent heat storage using solid-liquid phase change is more attractive when compared to other methods of thermal energy storage due to its large heat storage capacity at constant or nearly constant temperature corresponding to the phase change transition temperature of the PCM used (Dutt et al., 2004), therefore creating high storage density systems. Because the phase change requires some time to complete, therefore it becomes possible to use LHS to smooth temperature variations in the controlled system. The storage capacity through latent heat using a mass m of PCM is given by: where L is the latent heat of fusion per unit mass (J/kg). The use of phase change materials (PCMs) can be found in solar energy storage systems for water heating, green houses, heating and cooling of buildings, cooking and waste heat recovery systems (Atul et al., 2009). 1.2.2 Phase Change Materials Classification of phase change materials PCMs are at the center of LHESS technology. Upon storing heat, the PCM begins to melt when the phase change temperature is reached. The temperature then stays nearly constant until the melting process is finished. The heat stored during the phase change process of the material is called latent heat. PCM can store 5 to 14 times more heat per unit volume than sensible storage materials such as water, masonry, or rock (Hale et al., 6 1991). The PCM to be used in the design of thermal storage systems should possess desirable thermophysical and chemical properties. A large number of phase change materials (organic, inorganic and eutectic) are available in any required phase change temperature range. A classification of PCMs is given in Fig. 1.2. Figure 1.2: Classification of phase change materials Organic compounds Some examples of organic compounds are polyethylene glycol, fatty acid and paraffin. These compounds were discarded in the past because they were more costly than common salt hydrates, have somewhat lower heat storage capacity per unit volume and, possibly because of a bias against petroleum derivatives during the energy crises. It has now been realized that some of these compounds have strong advantages that outweigh these shortcomings. Some of the advantages are ability of congruently melting, self Eutectics Non- Paraffin Salt Hydrates Metallic Inorganic-Inorganic Inorganic-Organic Organic-Organic Inorganic Organic Paraffin Phase Change Materials 7 nucleating properties, non-corrosive behavior, chemically and thermally stable, but the disadvantages are low conductivity, low phase change enthalpy, and inflammability (Abhat et al., 1983). The normal paraffins of type CnH2n+2 are a family of saturated hydrocarbons with very similar properties. Paraffins between C5 and C15 are liquids, and the rest are waxy solids. Paraffin wax is the most used commercial organic heat storage PCM (Lane, 1983). It consists of straight chain hydrocarbons that have melting temperatures from 23 to 67°C (Abhat, 1983). Commercial grade paraffin wax is obtained from petroleum distillation and is not a pure substance, but a combination of different hydrocarbons. In general, the longer the average length of hydrocarbon chain, the higher the melting temperature and heat of fusion (Suwondo et al., 1994). Paraffins are easily available from many manufacturers and are usually more expensive than salt hydrates (Lane, 1983). Advantages Paraffin waxes show no tendency to segregate. They are also chemically stable although Lane (1983) reports slow oxidation when exposed to oxygen requiring closed containers. Stable properties after 1500 cycles in commercial grade paraffin wax have been reported (Liu et al., 2005). They also have no tendencies to super cool, so nucleating agents are not necessary (Lane, 1983). Paraffin waxes are safe and non-reactive (Hasnain, 1998). They are compatible with all metal containers and easily incorporated into heat storage systems. Care however should be taken when using plastic containers as paraffins have a tendency to infiltrate and soften some plastics (Lane, 1983). Disadvantages Paraffins have low thermal conductivity in their solid state. Paraffins have a volume change between the solid and liquid stages. This causes many problems in container design (Hasnain, 1998). Paraffins are flammable; this can be easily alleviated by a proper container (Suwondo et al., 1994). Lane (1983) also reports that paraffins can contract enough to pull away from the walls of the storage container greatly decreasing heat transfer to the system. Unlike salt hydrates, commercial paraffins generally do not have sharp well-defined melting points (Lane, 1998). 8 Inorganic compounds The popular inorganic compounds are the hydrated salts. They have a wide range of applications and some of the advantages associated with inorganic compounds are high volumetric latent heat storage capacity, often twice the capacity of organic compounds, high thermal conductivity and greater phase change enthalpy. The disadvantages are under cooling, phase separation, corrosion, phase segregation and lack of thermal stability (Abhat et al., 1983). Salt hydrates are some of the oldest and most studied heat storage PCMs (Lane, 1983). They consist of a salt and water, which combine in a crystalline matrix when the material solidifies. They can be used alone or in eutectic mixtures (Abhat, 1993). There are many different materials that have melting ranges from 15 to 117°C (Lane, 1983). Salt hydrates are the most important group of PCMs, which have been extensively studied for their use in LHESS. Three types of behavior of the melted salts can be identified: congruent, incongruent and semi congruent melting (Lane, 1998). Advantages Low cost and availability of salt hydrates makes them commercially attractive for heat storage applications. Two of the least expensive and most available salt hydrates are CaCl2∙6H2O and Na2SO4∙10H2O (Lane, 1983). Salt hydrates have a sharp melting point; this maximizes the efficiency of a heat storage system. They have a high thermal conductivity when compared with other heat storage PCMs, resulting in an increased heat transfer in and out of the storage unit. They have a high heat of fusion to decrease the needed size of the storage system. Salt hydrates have a lower volume change than other PCMs (Choi et al., 1992; Lane, 1999; Liu et al., 2005) Disadvantages These PCMS have a tendency to segregate, i.e., formation of other hydrates or dehydrated salts that tend to settle and reduce the active volume available for heat storage. Abhat reports a decrease in the heat of fusion of over 73% in Na2SO4∙10H2O after 1000 melt/freeze cycles (Abhat, 1993). This problem can be eliminated to a certain extent through the use of gelled or thickened mixtures (Lane, 1983) though this process 9 negatively impacts the heat storage characteristics of the mixture and the mixture still degrades with time (Abhat, 1993). Salt hydrates exhibit super cooling, resulting in a start of crystallization well below the freezing point of the PCM. This can be avoided using suitable nucleating materials to start crystal growth in the storage media (Choi et al., 1992). Another problem of salt hydrates is their tendency to cause corrosion in metal containers that are commonly used in thermal storage systems (Abhat, 1993), so compatibility of PCM and container should always be checked before use. Eutectic compounds A eutectic system is a mixture of chemical compounds or elements that has a single chemical composition and has a minimum melting composition corresponding to two or more components, each of which melts and freezes congruently forming a mixture of the component crystals during crystallization (Lane, 1999), leaving little opportunity for the components to separate. These compounds are not widely use; therefore, a more elaborate discussion on them is not necessary here. 1.2.3 Melting Phase Change Heat Transfer Melting heat transfer occurs in different geometries and many fields of engineering (Khillarka et al., 2000). For many years melting heat transfer has been a subject of intensive experimental, analytical and numerical studies, therefore the literature that concerned the melting phenomena is voluminous (Fukusaku et al., 1999). Phase change problems are also often called moving boundary problems or Stefan’s problem (Abhat et al., 1983) and their study presents one of the most exciting and challenging areas of current applied mathematical research (Lamberg, 2003). The existence of a moving boundary generally means that the problem does not admit a simple closed form analytical solution and accordingly much research has focused on approximate solution techniques (Khillarka et al., 2000). In general, phase change problems involve a transient, non-linear phenomenon with a moving liquid–solid 10 interface whose position is unknown a priori and also flow problems associated with thermal fluid. In addition, the two phases of PCM may have different properties and configuration of the LHESS unit may differ with the applications (Ghasemi et al., 1999). Conduction and natural convection are the dominant heat transfer mechanisms during the phase change process (Esam et al., 2004), although close contact melting plays an important role in certain configuration (Groulx et al., 2007). Classically, Stefan’s problem was first approached as pure conduction in semi-infinite medium (Lauardini, 1981) and later on natural convection has been considered during melting and solidification of PCMs. The different classes of solutions available for Stefan’s problem are analytical and numerical. Many approximate analytical techniques such as the heat balance integral, variation technique (Crank et al., 1995), isothermal migration, source and sink method (Buddhi et al., 1988) and periodic solution (Lazaridas, 1990) have been employed. A common drawback of these approximate techniques is the limitation to one dimensional analysis since they become very complex when applied to multidimensional problems. Numerical methods, both finite difference (Lazaridas, 1990) and finite element (Rolph et al., 1982) appear more powerful in solving the moving boundary problem. In general, a time variant mesh approach (Yimer et al., 1997) offers accuracy but is limited to simple problems and geometries. The fixed mesh approach (Yingqiu et al., 1999), in which the latent heat of fusion is usually absorbed into the material’s specific heat, or enthalpy, is much simpler and practical to use. Therefore, the numerical formulations widely used so far are the enthalpy and effective heat capacity methods (Lamberg, 2003). In the enthalpy method, the enthalpy, which is a function of temperature, is considered a dependent variable along with temperature. Thus the enthalpy based conduction equation is valid for both solid and liquid phases, as well as the solid–liquid interface. Therefore, there is no need to track the interface which makes this formulation attractive (Kim, 1992). In the effective heat capacity method, the heat capacity of the PCM during the phase change process is introduced. The effective heat capacity is directly proportional to the stored/released energy during phase change 11 and the specific heat. Therefore, with effective heat capacity, it is possible to describe non-isothermal phase change in the PCM. It is important to understand that governing heat transfer mechanisms may be different for different phase change processes (melting and solidification) and for different configurations, orientations of the system and fins. The following subsection will present a review of research work done in those various geometries. Rectangular Geometry Melting of phase change materials in rectangular enclosures has received considerable attention. It was demonstrated numerically that free convection plays a role during the melting process encountered in rectangular geometries (Wang et al., 1999; Khillarkar et al., 2000; Gong et al.,1998). The melting of an unfixed solid in a square cavity by using the fixed grid enthalpy method was studied (Ghasemi et al., 1998) and a numerical technique to describe the behavior of the flow at the boundary of the moving phase was also developed (Szimmat et al., 2002). Uros et al. (2003) worked on enhanced heat transfer in rectangular PCM thermal storage experimentally and reported that natural convection was present and increased heat transfer during melting and that fins could be used to effectively enhance melting and solidification. Asako et al. (1999) experimentally investigated the effect of the density change on PCMs melting rate in an unfixed rectangular enclosure under low gravity environment and Wang et al. (1998) looked at the melting process in a rectangular enclosures. Analytical and experimental studies found that convection played a role in the process of close contact melting of high Prantdl number substances (Groulx et al., 2007) and determined the behavior of ice during close contact melting on a sliding heated flat plate (Groulx et al., 2006). Mbaye et al. (2001) and Mohammed et al. (2003) analytically investigated the melting process of a solid PCM in a rectangular enclosure heated from a vertical side at a constant heat flux. 12 Spherical Geometry The melting process in a spherical geometry has attracted less attention until recently. Spherical capsules are now more widely used as storage of PCM in LHESS (Saitoh et al., 1998). Assis et al. (2007) numerically studied the melting process of a PCM in spherical geometry. They performed a detailed parametric investigation of melting in spherical shells of various diameters with a uniform temperature. Eames et al. (2002) performed experiments to determine the freezing and melting behavior of water in spherical enclosures of the types used in thermal storage systems and Moore et al. (1984) studied melting of an unconstrained PCM inside a sphere. The melting process of a PCM (n-octadecane) inside spherical capsules of various sizes and the melting in spherical capsule considering the effect of close contact and natural convection heat transfer was studied experimentally by Saitoh et al. (1993), and Saitoh et al. (1998) respectively. Ro et al. (1990) conducted an experiment looking at the melting of a PCM initially solid and at its fusion temperature contained in a spherical capsule. Also, Bahrami et al. (1987) analytically investigated melting in a spherical capsule and obtained a closed form distribution of temperature. Cylindrical Geometry Cylindrical enclosures are some of the most important geometric arrangements in view of their practical applications, such as in casting processes, and thermal storage systems (Benjamin et al. 2006). Numerical simulations based on either moving or fixed grids have been widely used in the study of solid/liquid phase change problems, facilitated by rapid increases in computational power (Brent et al., 1988; Simpson et al., 1998). Variants of the enthalpy method, such as the enthalpy-porosity method (Brent et al., 1988) and apparent heat capacity method (Simpson et al., 1998) have been employed. Numerical techniques based on moving grids have also been proposed (Zhang et al., 1996; Kim et al., 1992). A comparison of fixed and moving grids was performed by Viswanath et al. (1993) and later by Bertrand et al. (1999). They found that the appropriate choice in the solution method is often problem-dependent. Therefore, 13 advances in numerical modeling depend upon validation against rigorously controlled and well-documented experimental results. Experimental data for the melting of gallium and the solidification of tin (Wolff et al., 1988; Gau et al., 1986) have been employed by various investigators for validating predictions (Zhang et al., 1996; Kim et al., 1992). Benjamin et al. (2006) used photographic and digital image processing techniques and the enthalpy method to study the melting in a cylinder experimentally and numerically with the aim to provide benchmark experimental measurements of a PCM (n-ecosane) for validation of numerical codes. Also, experimental and analytical studies of the melting behavior of paraffin wax as a PCM in latent thermal energy storage using cylindrical capsule were performed; the enthalpy method was used to perform the analysis (Regin et al., 2006). Another work was carried out by Gong et al. (2008). They analyzed contact melting around a horizontal elliptical cylinder heat source and also shown that the melting process was mainly governed by the PCM melting temperature range and capsules radius. An experimental and numerical investigation of heat transfer during technical grade paraffin melting and solidification in a shell and tube latent thermal energy storage was also performed (Anica, 2005). Hendra et al. (2005) analytically and experimentally determined the thermal and melting heat transfer characteristics in a cylindrical latent heat storage system using PCM called Miikro (Indonesian traditional substance). Thermal Enhancer on LHESS Most PCMs have unsuitably low thermal conductivity, leading to slow charging and discharging rates, hence heat transfer enhancement techniques are required for most LHESS applications. Several studies have been conducted to study heat transfer enhancement techniques in PCMs using fins of different configurations (Ermis et al., 2007; Ismail et al., 2001; Choi et al., 1992; Agyenim et al., 2008; Sasaguchi et al., 1994; Zhang et al., 1996; Velraj et al., 1997). The majority of the heat enhancement techniques have been based on the application of fins embedded in the PCM. This is probably due to the simplicity, ease in fabrication and low cost of construction. The general observations 14 drawn from the various studies demonstrate that, irrespective of the PCM used, the heat transfer characteristics of the PCMs can be improved using all of the different enhancement techniques. Nayak et al. (2006) numerically studied the thermal performance of heat sinks with PCMs and thermal conductivity enhancers (porous matrix, plate-type and rod-type fins). They showed that in the case of a PCM combined with a porous thermal conductivity enhancer aluminum matrix, the thermal conductivity and melting rate were increased. Andrew et al. (2006) experimentally investigated the thermal conductivity enhancement of PCMs using a graphite matrix and they reported that graphite matrix impregnated with paraffin is a viable way of improving thermal conductivity. Choi and Kim (1992), Horbaniuc et al. (1999), Velraj et al. (1997) and Hamada et al. (2003), each employed different experimental setups, different container materials, and different PCMs to investigate the heat transfer enhancement characteristics of PCMs. In terms of performances of heat transfer enhancement techniques and systems used, the best enhancement technique as reported in the literature was that due to Velraj et al. (1997), where the effective thermal conductivity calculated employing paraffin with rings was ten times (2 W/m∙K) greater than the thermal conductivity of paraffin (0.2 W/m∙K). During the charging and discharging stages of the experiment, the thermal fluid in the experiment from Horbaniuc et al. (1999) passed over just a fraction of the whole setup rendering heat addition and removal inefficient. Results derived from any such setup cannot be considered to be optimal. Different researchers used different parameters to assess the heat transfer enhancement in the PCMs. Velraj et al. (1999) evaluated the enhancement of the heat transfer using the effective thermal conductivity taken from a two-dimensional enthalpy-temperature governing equation which assumed no variation of temperature and thermal conductivity in the axial direction. Results were presented graphically using temperature–time curve and as such limits the application tending to other applications. Horbaniuc et al. (1999) measured the performance of fins in terms of the interface freezing stage and the time taken for complete solidification to be achieved using parabolic and exponential approximations. 15 Hamada et al. (2003) used the effective thermal conductivity proposed by Fukai et al. (2002) to assess and compare results to the control system having no heat transfer enhancement. In the case of Choi and Kim (1992), the key parameters used to assess the heat transfer enhancement of the circular finned system were the ratio of overall heat transfer coefficient in the finned and the unfinned tube systems. They reported a ratio of overall heat transfer coefficient of 3.5 and 3.2 between the finned and the unfinned tube systems. A comparison was also made by deriving a relationship for the ratio of the total amount of heat recovered in the finned and unfinned tube systems correlated with dimensionless parameters of Fourier, Stefan and Reynolds numbers with a correlation coefficient of 0.994 and a standard deviation of 0.023 and 0.028, respectively. The above discussion illustrates the fact that there is no unified international or national standard methods (such as British Standards or EU standards) developed to test PCMs, making it difficult for comparison to be made to assess the suitability of PCMs to particular applications. Contradictions also exist in the thermophysical properties of PCMs provided, especially the latent heat of fusion, thermal conductivity and densities in solid and liquid states. This is again due to the absence of unified certification standards and procedures. However, based on the fact that finned systems have been reported as one of the most practical and easy to fabricate; fins were employed in this work to enhance the thermal transfer rate of the LHESS studied. 1.3 Problem Statement The goal of this research work is twofold: 1. To develop a numerical method of studying and solving phase change heat transfer problems encountered in a cylindrical latent heat energy storage device during charging (melting of the PCM) using the finite element method and the commercial software COMSOL Multiphysics. The main question to be answered is: How to account for the large amount of latent energy stored during melting of 16 a PCM, and the moving interface between the solid and liquid phase, in COMSOL using the finite element method? 2. As previously mentioned in the literature review, PCMs have low thermal conductivity. To remedy the situation, the heat transfer characteristics of a LHESS can be enhanced by properly designing such a system. This study looks at the effect of adding fins on the charging (melting) rate of the PCM, and subsequently investigates the effect of the thermal fluid velocity on the melting rate of the PCMs. Due to the inherent non-linearity of the energy balance at the solid-liquid interface governing the interface behavior during melting, few analytical solutions of interest are available, giving rise to the development of a great number of numerical algorithms based on the finite differences, finite elements, and more recently, boundary elements methods (Atul et al., 2009). The algorithm usually account for the melting interface movement in two ways: the front tracking or fixed domain methods. The main disadvantage of front tracking methods is that they require some degree of irregularity in the treatment of the moving boundary and its evolution. In contrast, fixed domain methods, which are based on weak formulations of the energy equation, can handle complex topological evolutions of the interface (Zivkovic et al., 2001). Furthermore, the association of fixed domain techniques with finite elements is quite natural due to their common search for the modeling of complex geometries. In this study, COMSOL Multiphysics functions are used to account for the discontinuity at the interface. The treatment of the energy balance at the melting interface, using a fixed domain method, call the effective heat capacity (Lamberg, 2003) or enthalpy methods (Date, 1992). In solving the phase change problem numerically with COMSOL Multiphysics, the heat capacity method is studied and adopted; the method is also validated analytically. Finally, the present study focuses on accounting for the total energy stored in a LHESS, for that reason, conduction in the PCM is the only mode of heat transfer considered; making the present problem akin to solving a multidimensional Stefan’s problem (moving boundary problem). 17 1.4 Scope of the Thesis Chapter 2 presents the LHESS geometry analyzed in this work as well as the materials used. The physical processes encountered during the charging process and the related governing equations with the initial and boundary conditions are also presented. It is shown how the physical equations and boundary conditions are treated numerically using finite element analysis in COMSOL Multiphysics. Chapter 3 presents the validations of the three processes present in the LHESS study: fluid flow, heat transfer and phase change. The numerical results are validated using developed analytical solutions. In order to properly simulate the physical processes involved, certain steps were taken in performing the numerical simulations which involve a mesh convergence study as well as selecting the appropriate time steps. These considerations are also presented. The LHESS numerical results are summarized in chapter 4. The effects of the heat transfer fluid (HTF) inlet velocity, and of fins, on the phase change process are discussed. The total energy stored in the LHESS after 12 hours of charging is also calculated for various fin configurations and the LHESS efficiency introduced. In chapter 5, concluding remarks are presented and recommendations for future studies on this problem are given. Appendices A, B and C present the temperature plots obtained as a function of the thermal fluid velocity for various number of fins, number of fins with thermal fluid inlet velocity of 0.05 m/s and number of fins with thermal fluid inlet velocity of 0.5 m/s respectively. 18 CHAPTER 2: MODEL GEOMETRY AND MATHEMATICAL FORMULATION This chapter presents the geometry of the LHESS studied in the numerical simulations and the materials employed in the LHESS. The arrangement of the annular fins in the storage device is also presented. The physical processes encountered in the thermal fluid study of the LHESS, as well as the mathematical formalism that need to be solved, i.e., equations and boundary conditions, are then presented. Information’s on how to treat the mathematical formalism numerically in COMSOL Multiphysics are also shown. 2.1 LHESS Geometry and Materials The geometry of the cylindrical LHESS used in this study is presented in Fig. 2.1. It is a cylinder having a copper pipe running through its center. A thermal fluid flows in the pipe while the rest of the cylinder is filled with PCM. The dimensions of the LHESS are given in Table 2.1. The list of materials used, as well as their thermophysical properties is presented in Table 2.2. Notice the low thermal conductivity of paraffin wax (0.21 W/m∙K) compared to that of copper (400 W/m∙K), as well as the high latent heat of fusion and heat capacity of the PCM, which are all critical in the charging process of the LHESS. The radial fins are evenly spaced. However, in order to allow for uniform melting of the PCM between the compartments created by the fins, the first and last compartments, bound between the first and last fin, and cylinder end walls, are half the size of compartments formed by two adjacent fins. This must be done because less heat is transferred into the end compartments formed by the two end fins and the cylinder wall, compared to the heat transferred in compartments found between two adjacent fins. 19 (a) (b) Figure 2.1: (a) Geometry of cylindrical LHESS; (b) Some examples of fin configurations in the cylindrical LHESS (0 to 3 fin arrangement are shown) Table 2.1: Materials and Dimensions of the LHESS Geometry Material Length (mm) Inner radius (mm) Outer radius(mm) Thickness (mm) Shell Insulation 1000 303 313 10 Tube Copper 1400 27 30 3 Fin Copper - 30 303 5 20 Table 2.2: Thermophysical Properties of Materials Materials Density(kg/m 3 ) Thermal Conductivity (W/m∙K) Heat Capacity (kJ/kg∙K) Latent Heat of Fusion (kJ/kg) Copper 8700 400 0.385 - PCM (Paraffin wax) 750 0.21 2.4 174 Insulation 1150 0.26 1.7 - 2.1.1 Fins Selection The radial fin selected is 5 mm thick and extend the full radius of the LHESS. The thickness is to achieve reasonable strength, in order for the fins to support the weight of the PCM place in between them. However, a study of the effect of the fin thickness was not part of this study. Emphasis was laid on the effect of addition of fins on the LHESS as stated in the research objective. 2.2 Physical Processes The four physical processes involved in the thermal fluid behavior of the LHESS are: fluid flow, heat transfer by conduction and convection, and phase change heat transfer. This section presents the general governing equations and the corresponding initial and boundary conditions needed to solve for each physical processes. 2.2.1 Fluid Flow Process Considering incompressible axisymmetric flow, the fluid dynamics governing equations are the Navier-Stokes and the continuity equations. These equations have to be solved simultaneously, in order to account for the fluid flow process encountered within the working fluid in the LHESS. The axisymmetric continuity equation for an incompressible flow is given as (Deborah et al., 2005): 21 While the incompressible axisymmetric Navier-Stokes equations are (Deborah et al., 2005): These equations must be solved using the following initial and boundary conditions: Where u0 is the initial inlet fluid velocity and is the inner radius of the copper pipe used. Relative pressures are also calculated as a function of an absolute value of pressure defined at a specific point in the flow field. 2.2.2 Heat Transfer Process Heat Transfer by Conduction Heat transfer in the PCM is by conduction only, since the effect of convection in the melted PCM is assumed to be negligible, due to the small volume of PCM between the fins of the system for higher number of fins configurations and also the non movement of the PCM from one compartment to the other. The general axisymmetric energy equation in cylindrical coordinate is (Deborah et al., 2005): 22 and needs to be solved using the following initial condition: As for boundary conditions, every outside surface of the system is insulated and must respect the following condition: , where represents a flow normal to any outside surface. And, the heat flux through the inner surface (copper pipe) is determined by continuity from the heat flux by convection calculated in the fluid flow. Heat Transfer by convection in PCM During the melting process, heat is transferred from the wall to the phase change material first by conduction, and then by natural convection. This is because, as the solid region moves away from the heat transfer surface and the thickness of the liquid region increases, natural convection effects become more pronounced. The influence of natural convection on the location of the melt front has been investigated by Lamberg et al. (2003). The melting in a semi-infinite PCM storage, with an internal fin, was studied analytically. The analytical model used the well-known Newmann solution, which assumes that the heat is transferred only by conduction, neglecting natural convection (Jegadheeswaran et al., 2009; Eva et al., 2007). It was found that although the Newmann solution is exact, it underestimates the location of the melt front. Again, Lamberg et al. (2004) conducted a numerical study on the melting of a PCM inside a rectangular enclosure with and without the effect of natural convection, the dimension of the rectangular storage was 96 mm high (equivalent to the height of each compartments in the 10 fin configuration) and 20 mm wide. The results were compared to experimental results. It was observed from the results that when the natural convection effect was ignored, the PCM took longer to reach the maximum temperature, which shows the limitation of ignoring natural convection larger PCM compartment. When it comes to natural convection, the Rayleigh number is the dimensionless number guiding the physical process. The critical Rayleigh’s number for the onset of structured cell-like convection in a horizontal rectangular cavity heated from below is: 23 Table 2.3: The PCM Fluid Properties Material Thermal diffusivity (m 2 /s) Kinematic viscosity (m 2 /s) Thermal expansion coefficient (K -1 ) Paraffin wax 1.16 × 10 -7 6.61 × 10 -6 1.1 × 10 -4 When the Rayleigh number is below the critical value of 1708, heat transfer is by conduction. Convection in the form of fluid motion consisting of regularly spaced roll cells in the rectangular cavity appears when the Rayleigh number exceeds 1708. The geometry studied in this work can be represented by radial fins in a cylindrical geometry, since most of the heat transfer is by the fins, with the cold surface being the solid-liquid melting interface which is non-uniform and constantly varying with time. This geometry does not conform directly to the applied correlation. For this reason, it is doubtful that structured convection cells would be present in the system. Unstructured natural convection, the type that is expected to be present in the studied system, appears for Rayleigh numbers above which correspond, in this case, to a characteristic length (vertical half-distance between two fins) of 2 cm (Eq. (2.11)). This characteristic length is of the order of the half-distance between two fins for systems having more than 12 fins (the half-distance between fins keeps on decreasing with addition of fins). with Tw = 350 K and T∞ = 316 K. Based on this order of magnitude calculation, natural convection effect reduces with an increase in the number of fins to the system, since as the number of fins added reduces the size of the PCM compartment. Convection would play a greater role in the configurations having less than 12 fins, even though the amount of melted PCM in these systems has a thickness of the order of 2 cm in some cases. For configurations having more than 12 fins, convection would play a very small part during melting, the effects of convection in these cases being more pronounced when the entire PCM is melted. 24 Heat Transfer by convection in Thermal Fluid Convection is the prevailing heat transfer process within the thermal fluid to the wall of the pipe. The energy equation that needs to be solved in that case is given by (Deborah et al., 2005): Where represents the heat generated by viscous dissipation, which has an insignificant effect in the present problem. This equation, coupled to the Navier-Stokes and continuity equations, must be solved using the following initial and boundary conditions: The heat flux calculated at the wall of the copper pipe is directly related to the heat flux used to solve the heat conduction equation in the rest of the LHESS. 2.2.3 Phase change heat transfer During melting of the PCM, energy is stored in the LHESS. Therefore it is important to account for the phase change process taking place within the PCM. The general equation for the phase change heat transfer process at the melting interface is given by (Naterer, 2003): Where refers to the position of the solid-liquid interface. During melting, at this solid- liquid interface, the energy balance requires that the difference between the heat flux in the solid and liquid phases be balanced by the latent heat absorbed during phase change. In the LHESS studied, the entire PCM is solid initially, and will melt when its temperature reached the melting point. 25 2.3 Numerical Method The finite element method (FEM), or finite element analysis (FEA), is one of the most popular numerical computing techniques used in the field of engineering. This section describes how a finite element software, COMSOL Multiphysics, is use to simulate the systems of linear and non-linear time dependent partial differential equations (PDEs) encountered in the physical description of the thermal fluid process found in the studied LHESS. 2.3.1 Numerical Modeling After creating the geometry, in this case a 2D axisymmetric model (a simple 2D model is used in the validation section), the initial and boundary conditions are applied and then a mesh is created by discretizing the computational domain into triangular mesh elements, i.e., partitioning the geometry into small units of simple shapes called mesh elements. In order to properly simulate the finite element model, care is taken in selecting the time step. Finally the numerical simulation is performed using COMSOL Multiphysics. The fluid flow and the heat transfer (conduction and convection) equations are built directly into COMSOL Multiphysics, therefore the software solves these systems of non-linear PDE simultaneously by applying the initial and boundary conditions in each case. The following subsections present the applied initial and boundary conditions for the fluid flow and heat transfer processes simulated in this study, as they are defined in COMSOL Multiphysics. Fluid Flow applied boundary conditions Inlet boundary condition: This boundary condition defines the normal inflow velocity of the thermal fluid in the LHESS. This is expressed as (COMSOL Multiphysics User’s Guide, 2008): Where is the boundary normal which points out of the domain. 26 Wall boundary condition: This boundary condition describes the existence of a solid wall, hence a no slip condition resulting in a thermal fluid velocity of zero at the wall. It is expressed as (COMSOL Multiphysics User’s Guide, 2008): Outflow boundary condition: This boundary condition is used to define the thermal fluid behavior exiting the LHESS. Pressure, no Viscous Stress is the applied boundary condition in this case. This boundary condition prescribes vanishing viscous stress along with a Dirichlet condition on the pressure. It is a numerically stable boundary condition that admits total control of the pressure level at the whole boundary. It is expressed mathematically as (COMSOL Multiphysics User’s Guide, 2008): Symmetry boundary condition: This boundary condition is applied by taking advantage of the symmetry of the system. It results in the use of a 2D axial symmetry model instead of a complete 3D geometry. This saves memory space and reduces computational time. Heat Transfer applied boundary conditions Inlet boundary condition: This boundary condition defines the thermal fluid inlet temperature. The inlet boundary condition is expressed as (COMSOL Multiphysics User’s Guide, 2008): Continuity boundary condition: The interior boundaries of the LHESS are defined by these boundary conditions. It allows the heat flux across the internal boundaries of the system in the normal direction to be continuous. Also the temperature is naturally continuous along the internal boundaries due to the continuity of the finite element field 27 linking neighboring elements. The Continuity boundary condition is identical to the condition that applies between any two neighboring elements in the mesh which allows for continuous heat flow between them. Mathematically it is given as (COMSOL Multiphysics User’s Guide, 2008): Where q1 and q2 are the calculated heat flux on two neighboring elements and Eq. (2.22) defines these heat fluxes. Insulation boundary condition: In order to prevent heat loss from the LHESS, this boundary condition is used where the domain of the LHESS would be insulated in a real- life situation. Intuitively the insulation boundary equation states that the thermal gradient across the boundary must be zero. The boundary equation is expressed as (COMSOL Multiphysics User’s Guide, 2008): Phase Change Modeling The phase change phenomenon is modeled separately due to the non-linear nature of the problem. A wide range of different numerical methods for solving PCM problems exist, as were presented in chapter 1. The major problem associated with the numerical analysis of the phase change process is the need to account for the melting interface position and the large amount of energy needed to melt the PCM. This problem is dealt with by using the effective heat capacity method. During the melting process, the effective heat capacity of the PCM is used. A discontinuous function defined by Eq. (2.24), accounts for the total amount of energy stored in the PCM. The effective heat capacity, Cp,eff of the PCM is calculated using the following equation (Lamberg, 2003): 28 Where L is the latent heat of fusion (J/kg), T1 is the onset temperature of phase transition, T2 is the temperature when the phase transition ends, Cps and Cpl are the respective solid and liquid phase heat capacities. In the present study, Cps and Cpl have the same value, therefore, the variable Cp will be used for both. For the PCM used in this study, paraffin wax, the heat capacity is modified and written in the following way (Lamberg, 2003): Discontinuous functions are created in the COMSOL Multiphysics in order to solve the solid-liquid phase change problem numerically. The Cp of paraffin wax, presented in Eq. (2.25), is incorporated by using a continuous step function defined as (COMSOL Multiphysics User’s Guide, 2008): Where flc2hs is a smoothed Heaviside function with a continuous second derivative without overshoot. Figure 2.2 shows the variation of the heat capacity with temperature of the paraffin wax used in the simulation when Eq. (2.26) is applied numerically (Eq. (2.27)) (COMSOL Multiphysics User’s Guide, 2008). This specific heat presented clearly shows the modified specific heat capacity which accounts for the high specific heat encountered within the melting temperature range, during phase change. This accounts for the large amount of energy (latent heat) stored during melting of the PCM. 29 Figure 2.2: Modified heat capacity of paraffin wax as a function of temperature. 0 10 20 30 40 50 60 70 290 300 310 320 330 340 350 H e a t C a p a ci ty o f P C M ( k J/ k g K ) Temperature (K) 30 CHAPTER 3: VALIDATION The description of the governing equations needed for the physical processes encountered in the LHESS has been given in the previous chapter. In this chapter, the analytical validation of the specific numerical results obtained for the thermophysical processes encountered in the LHESS will be discussed. Section 3.1 describes the validation of the fluid flow simulations performed on the heat transfer fluid in the pipe of the LHESS; section 3.2 presents the validation of the heat transfer simulations performed within the LHESS; while the validation of the phase change process that takes place within the PCM is presented in section 3.3. Mesh convergence studies are carried out for every numerical simulation in order to guarantee accurate numerical solutions to every thermophysical process studied. 3.1 Fluid Flow Validation 3.1.1 Analytical Solution The fluid flow problem solved to validate the capacity of COMSOL Multiphysics to accurately address such force flow problem is the simple 2D Poisseuille flow presented in Fig. 3.1. Figure 3.1: Schematic representation of Flow between two plates x y u (y) 31 The incompressible flow is forced between two horizontal, infinite parallel plates by a constant pressure gradient forcing the fluid to move in the x-direction, parallel to the plates. To solve this problem in Cartesian coordinates, the following continuity equation must be solved (Deborah et al., 2005): along with the steady-state Navier-Stokes equations for the x-component of momentum (effect of gravity neglected) (Deborah et al., 2005): The boundary conditions needed are as follow: Using the previous boundary conditions, Eq. (3.2) can be rewritten, by accounting for Eq. (3.1), as follow: This can be integrated twice to profile an expression for the velocity profile between the plates: Applying the no-slip boundary conditions (Eq. (3.3)), a parabolic velocity distribution is obtained: This can be restated in a different form using the average velocity calculated as: 32 resulting in the following velocity profile (Deborah et al., 2005): 3.1.2 Numerical Analysis and Validation The 2D Poisseuille flow, for which an analytical solution was presented in section 3.1, was solved numerically. Figure 3.2 presents the 2D numerical model geometry used to solve for this flow. A thermal fluid (water) flows through a 0.1 m thick region form by two fixed parallel copper plates of equal thickness and length of 0.3 m and 50 m respectively. Such a long system was used in order to deal with the dynamic entrance region and obtain a flow as close as possible to a fully developed flow. The applied initial and boundary conditions are presented in Table 3.1. The initial system condition is the condition of the entire system before the boundary conditions are applied. Figure 3.2: Numerical model geometry and finite element mesh used to solve for a flow of water between two parallel plates. No slip condition Thermal Insulation 33 Table 3.1: Initial and boundary conditions for the numerical fluid flow model. Initial system conditions Inlet conditions (water) Outlet conditions (water) Channel wall conditions Fluid 0.05 m/s 0.05 m/s No viscous stress P0 = 0 No slip Simulations were performed for four different element sizes (0.12 m, 0.08 m, 0.04 m and 0.01 m), and for a simulated time of 28 hours. The element size, number of elements, and degree of freedom (DOF) used for the numerical simulations, as well as the total COMSOL Multiphysics simulation time resulting from each simulation, are presented in Table 3.2. Table 3.2: Element sizes, number of elements, DOF and simulation time for the numerical fluid flow model. Element size (m) Number of elements Number of degree of freedom solved for Simulation time (seconds) 0.12 6 266 19 902 27 0.08 17 966 51 668 45 0.04 55 774 154 191 126 0.01 870 142 2 360 523 954 Figure 3.3 shows the position of the eight sections where velocity profiles obtained from the simulations were presented and analyzed. The positions of the sections are measured from the fluid inlet position (indicated by the arrow on Fig. 3.3), which corresponds to 9 m, 15 m, 25 m, 35 m, 43 m, 45 m, 47 m and 49 m and represented by sections a-a, b-b, c-c, d-d, e-e, f-f, g-g and h-h respectively. 34 Figure 3.3: Position of the presented and analyzed velocity profiles. Figure 3.4: Simulated fluid velocity profile, for various element sizes, taken on section h-h. Simulated time: 28 hours. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 1 1.02 1.04 1.06 1.08 1.1 F lu id v e lo ci ci ty ( m /s ) Distance (m) 0.12m 0.08m 0.05m 0.01m a b c d e f g h a b c d e f g h 35 Figure 3.4 presents the numerical velocity profiles obtained as a function of the element size used to model the geometry of the system. The results were obtained after 28 hours of simulated time and are taken at section h-h (49 m length) as shown in Fig. 3.3. The figure clearly shows that as the element size is reduced from 0.12 m to 0.05 m, the fluid velocity calculated increases. Further reduction of the element size to 0.01 m does not change the velocity significantly. Further reduction of the element size below 0.01 m resulted in system failure. The system ran out of memory during mesh processing. The mesh convergence result is presented in Fig. 3.5. The maximum fluid velocity is plotted against the element size. Reducing the element size below 0.05 m barely changes the maximum velocity obtained numerically, which is 0.071 m/s. Therefore, for maximum accuracy, element size of 0.01 m (870 142 total elements) was adopted for the numerical simulation and resulting validation. Figure 3.5: Maximum fluid velocity as a function of the element size used for the numerical fluid flow model at length 49 m. Simulated Time: 28 hours. 0.06 0.062 0.064 0.066 0.068 0.07 0.072 0.01 0.03 0.05 0.07 0.09 0.11 M a x im u m f lu id v e lo ci ty ( m /s ) Element size (m) 36 Based on the element size selected (0.01 m), a simulation was performed and the velocity profiles calculated at the different sections presented in Fig. 3.3 are shown Fig. 3.6. The development of the flow along the length of the system can be clearly observed from the presented velocity profiles, the flow becoming fully developed around the 43 m mark, since no significant changes are observed in the velocity profile after that point. Figure 3.6: Simulated fluid velocity profiles at different section. Simulated time: 28 hours. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 1 1.02 1.04 1.06 1.08 1.1 F lu id v e lo ci ty ( m /s ) Position(m) 9m 15m 25m 35m 43m 45m 47m 49m 37 The analytical velocity profile calculated from Eq. (3.9) is compared to the simulated fully develop numerical velocity profile in Fig. 3.7. Both analytical and numerical velocity profiles are obtained based on the system thickness (2 and thermal fluid inlet velocity of 0.05 m/s. It can be seen from Fig. 3.7 that both velocity profiles tend to agree, but a difference of 5% exists between the maximum velocity calculated numerically (0.071 m/s) and obtained analytically (0.075 m/s). This difference can be explained by the system limitation in reducing the element size used below 0.01 m in the numerical model; further mesh refinement would have resulted in closer results. Figure 3.7: Comparison of the analytical and numerical velocity profiles. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 1 1.02 1.04 1.06 1.08 1.1 F lu id i n le t ve lo ci ty ( m /s ) Position(m) Analytical Numerical 38 3.2 Heat Transfer Validation The validation of the heat transfer process focuses mainly on the convection heat transfer from the thermal fluid to the pipe wall. Using well known analytical methods and correlations, the heat transfer coefficient in the system will be calculated and the results compared with the numerical solution. 3.2.1 Analytical Solution Pure analytical solutions are almost impossible to obtain, even for the simplest geometries, for turbulent flows. Instead a combination of simplifying assumptions and experimental correlation were used to develop an analytical solution for the turbulent thermal flow studied in this section. Turbulent convection in pipe flow and turbulent flow over a plate were the approaches used to calculate the heat transfer coefficient. Turbulent Convection in Pipe Approach Gnielinski’s correlation is used for convective heat transfer (Deborah et al., 2005): which is valid for and . The thermal fluid properties are calculated at the mean temperature ( ). The flow Reynolds number must be calculated: as well as the turbulent flow friction factor (f) for circular tubes, calculated with Haaland’s equation (Deborah et al., 2005): 39 where is the copper pipe relative roughness and is the pipe absolute roughness taken as 1.5×10 -6 m for copper. The Prandtl number is defined as: Once the NuD is obtained from Gnielinski’s correlation (Eq. (3.10)), the convection coefficient (h) is obtained from (Incropera et al., 2007): Turbulent Flow over a Plate with uniform surface heat flux Approach This approach involves using correlations obtained from solving the boundary layer equations for a flow over a plate with a uniform surface heat flux imposed. In this case, all the fluid properties are evaluated at the mean boundary layer temperature ( ), also called the film temperature (Incropera et al., 2007): where Tw and T∞ are the wall temperature and the free-stream temperature respectively. Assuming a uniform surface heat flux imposed at the plate, the turbulent flow correlation stated in Eq. (3.16) is used to evaluate the local Nusselt number ( ) and the resulting convection coefficient (hx), where x is taken as the distance from the start of the flat plate (Incropera et al., 2007): In this case, the local Reynolds number is calculated as: 40 3.2.2 Numerical Analysis Figure 3.8 presents the 2D axisymmetric LHESS system studied, which is similar in dimensions and materials used as the one presented in chapter 2. The numerical simulation is performed using the applied initial and boundary conditions shown in Tables 3.3. Table 3.3: Initial and boundary conditions for the heat transfer study. Initial system conditions Inlet conditions (water) Outlet conditions (water) Channel wall conditions Outer surfaces of the plates Fluid 0.05 m/s 0.05 m/s No viscous stress, P0 = 0 No slip - Thermal 293K 350 K Convective heat flux Continuity Thermally insulated Reliable numerical solutions are obtained by performing a mesh convergence study. This is done by changing the size of elements used (mesh refinement) and running the simulation again until the numerical results converge to a stable value which is mesh independent. The numerical simulations were performed for five different element sizes and for a simulated time of 20 hours. Table 3.4 lists the five element sizes as well as the number of elements used, DOF and the time COMSOL Multiphysics took to run each simulation. In order to compare the calculated numerical convection coefficient (h) in the thermal fluid to the analytical convection coefficients, three equally spaced cross-sections in the pipe are chosen. They all fall within the 1 m length of the LHESS and are chosen as 0.45 m, 0.7 m and 0.95 m measured from the fluid inlet of the pipe (0.2 m before the beginning of the LHESS), which corresponds to section B-B, C-C and D-D respectively as shown in Fig. 3.9. The mesh convergence study is performed at section A-A shown in the figure. Figure 3.9 also shows the temperature obtained after 20 hours of simulated time. 41 Figure 3.8: Numerical model geometry and finite element mesh used to solve convection heat transfer in a pipe. Figure 3.10 presents the temperature profile along section A-A obtained numerically with the different element size used. The temperature is mesh dependent for element sizes of 0.018 m, 0.016 m and 0.012 m, but becomes mesh independent for element sizes of 0.008 m and 0.005 m. Further reduction of the element size below 0.005 m does not change the numerical results. Copper PCM PCM PCM PCM PCM PCM PCM PCM PCM PCM PCM PCM PCM PCM Insulation 42 Table 3.4: Element sizes, number of elements, DOF and simulation time for the numerical heat transfer model. Element size (m) Number of elements Number of degree of freedom solved for Simulation time (seconds) 0.018 18 985 50 934 8 456 0.016 19 257 52 530 9 035 0.012 19 843 56 649 9 412 0.008 20 759 60 181 9 864 0.005 36 782 95 929 12 658 Figure 3.9: Simulated temperature plot for the numerical model. Simulated time: 20 hours. D D Z C C Y X B B A A 43 The result of the mesh convergence study is shown in Fig. 3.11. The result indicates that a constant minimum temperature in the system is obtained for element sizes of 0.008 m and 0.005 m, making the solution mesh independent for these element sizes. To limit the solution times, an element size of 0.008 m was used in the numerical simulations. Figure 3.10: Temperature profile, for various element sizes, at section A-A. Simulated Time: 20 hours. 325 330 335 340 345 350 0.00 0.05 0.10 0.15 0.20 0.25 0.30 T e m p e ra tu re (K ) Distance from the system centerline (m) .018m 0.016m 0.012m 0.008m 0.005m 44 Figure 3.11: Minimum temperature as a function of the element size, at section A-A. Simulated time: 20 hours. The thermal fluid velocity and temperature profiles at the three sections for 20 hours of simulated times are presented in Figs. 3.12 and 3.13 respectively. It can be seen from Fig. 3.12 that the thermal fluid velocity and the velocity boundary layer thickness both increase along the pipe length since the fluid flow is still developing. Fig. 3.13 shows that the temperature of the fluid within the thermal boundary layer decreases along the pipe length but the thermal boundary layer thickness increases with length in the fluid flow direction. This is expected since the total energy carried by the thermal fluid decreases constantly along the length of the pipe since more and more energy is transferred to the rest of the system. 328 329 330 331 332 333 334 335 0.005 0.007 0.009 0.011 0.013 0.015 0.017 M in im u m T e m p e ra tu re ( K ) Element size (m) 45 Figure 3.12: Simulated thermal fluid velocity profiles at sections B-B, C-C and D-D. Simulated time: 20 hours Fluid Reference Temperature Figure 3.14 presents the schematic thermal fluid temperature profile obtained at sections B-B, C-C and D-D. In order to calculate the mean temperature at each section for various times, the profile is approximated as two straight lines. The obtained fluid reference temperatures are used to evaluate the thermophysical properties of the fluid. 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0 0.005 0.01 0.015 0.02 0.025 T h e rm a l fl u id v e lo ci ty ( m /s ) Distance from the system centerline (m) Velocity-0.95m Velocity-0.7m Velocity-0.45m 46 Figure 3.13: Simulated thermal fluid temperature profiles at sections B-B, C-C and D-D. Simulated Time: 20 hours. Figure 3.14: Schematic thermal fluid temperature profile along sections B-B, C-C and D-D. 325 330 335 340 345 350 355 0 0.005 0.01 0.015 0.02 0.025 T e m p e ra tu re ( K ) Distance from the system centerline (m) Temperature-0.95m Temperature-0.7m Temperature-0.45m X X Y r r1 r2 B A T Tmax Tmin Tmax=350K 47 Under this approximation, the mean temperature is calculated using the following equation: With: s is the slope of the angled straight line (Fig. 3.14) calculated as follow: and b is calculated knowing that: Integrating Eq. (3.19) yields the mean temperature used as the fluid reference temperature: The thermal fluid reference temperatures at the three aforementioned sections are calculated using Eq. (3.23) for simulated time of 5, 10, 15 and 20 hours; the results are presented in Table 3.5. Correspondingly, the pipe wall temperatures and heat flux obtained at those three sections (positions x, y and z on Fig. (3.9)) for the same simulated times are presented in Tables 3.6 and 3.7 respectively. Then, knowing , Tw and , and using Newton’s cooling law (Incropera et al., 2007): The convection coefficient, stemming from the numerical simulations, can be calculated: 48 Table 3.5: Fluid reference temperatures cross section B-B, C-C and D-D. Time (hours) Thermal fluid Reference Temperatures, (K) @ 0.45m @ 0.70m @ 0.95m 5 347 346 344 10 347 347 345 15 348 347 346 20 349 348 346 Table 3.6: Pipe wall temperatures at point x, y and z. Time (hours) Pipe Wall Temperatures K @ 0.45m @ 0.70m @ 0.95m 5 328 323 327 10 332 326 325 15 335 329 329 20 338 332 330 Table 3.7: Heat fluxes at points x, y and z. Time (hours) Heat Flux at the Wall ) W/m2 @ 0.45m @ 0.70m @ 0.95m 5 6110 6129 4490 10 5309 5612 4215 15 4649 5107 3945 20 3903 4753 3772 Table 3.8 and Fig. 3.15 present the numerical convection coefficient calculated, at different time and at the three sections, using Eq. (3.25). It can be seen in Fig. 3.15 that at any given time the heat transfer coefficient decreases in the flow direction due to the constant increase of the thermal boundary layer. However, the convection coefficient increases with time, irrespective of the section, with an average of 6% for every 5 hours 49 interval. A maximum numerical convection coefficient of 378 W/m 2∙K is obtained at section A-A (0.45 m) after 20 hours of simulated time. Table 3.8: Numerical convection coefficients (h) Time (hours) Values of convection coefficients (h) W/m 2 K @ 0.45m @ 0.70m @ 0.95m 5 322 263 201 10 335 274 212 15 345 282 219 20 378 301 234 Figure 3.15: Numerical convection coefficients as a function of time at sections B-B, C-C and D-D. 200 220 240 260 280 300 320 340 360 380 400 5 7 9 11 13 15 17 19 C o n ve ct io n c 0 e ff ic ie n t (W /m 2 K ) Time (hours) 0.45m 0.7m 0.95m 50 3.2.3 Validation Turbulent convection in a pipe approach Equations (3.10) to (3.14) were used to calculate the Reynolds numbers, friction factor, Prandtl number, Nusselt numbers and the convection coefficients needed to validate the numerical results using this approach. The calculated values, for 5, 10, 15 and 20 hours, are presented in Tables 3.9 to 3.13. Table 3.9: Calculated parameters using turbulent convection in pipe approach. Simulated Time: 5 hours. Pipe Length(m) Temperatures (K) f 0.45 328 0.0370 5358 3.22 32 0.70 323 0.0380 4882 3.58 30 0.95 318 0.0390 4498 3.93 28 Table 3.10: Calculated parameters using turbulent convection in pipe approach. Simulated Time: 10 hours Pipe Length(m) Temperatures (K) f 0.45 332 0.0365 5576 3.08 32 0.7 326 0.0375 5100 3.41 31 0.95 321 0.0384 4715 3.73 29 Table 3.11: Calculated parameters using turbulent convection in pipe approach. Simulated Time: 15 hours. Pipe Length(m) Temperatures (K) f 0.45 334 0.0361 5799 2.94 33 0.7 329 0.0370 5358 3.22 32 0.95 324 0.0378 4971 3.50 30 51 Table 3.12: Calculated parameters using turbulent convection in pipe approach. Simulated Time: 20 hours. Pipe Length(m) Temperatures (K) f 0.45 338 0.0356 6067 2.79 34 0.7 331 0.0365 5575 3.08 32 0.95 327 0.0373 5187 3..34 31 Table 3.13: Analytical convection coefficients using turbulent convection in pipe approach. Time (hours) Values of convection coefficients (h) W/m 2 K @ 0.45m @ 0.70m @ 0.95m 5 384 359 337 10 393 370 349 15 405 379 363 20 417 388 376 Figure 3.16 shows the comparison between the numerical and analytical (turbulent convection in pipe approach) convection coefficients for the same sections and times. Both numerical and analytical results clearly show that the heat transfer coefficient increases gradually with time for all sections. The percentage differences between the numerical and analytical convection coefficient at the 0.95 m, 0.7 m and 0.45 m sections are 40%, 26% and 16% respectively at any given time. Turbulent Flow over Plate with uniform surface heat flux approach Equations (3.15) to (3.18) were used to calculate the film temperatures, Reynolds numbers, Prandtl number, Nusselt numbers and the analytical convection coefficients for the sections at 0.45 m, 0.7 m and 0.95 m and for 5, 10, 15 and 20 hours. The results are presented in Tables 3.14 to 3.18. 52 Table 3.14: Calculated analytical parameters using turbulent flow over plate with uniform surface heat flux approach. Simulated Time: 5 hours. Pipe Length(m) Temperatures (K) 0.45 339 51353 2.75 242 0.7 336 76956 2.86 340 0.95 336 103543 2.89 433 Figure 3.16: Analytical (turbulent convection in pipe approach) and numerical convection coefficients at different sections as a function of time. 200 220 240 260 280 300 320 340 360 380 400 420 440 5 7 9 11 13 15 17 19 C o n ve ct io n c o e ff ic ie n t (W /m 2 K ) Time (hours) Num-0.45m Num-0.70m Num-0.95m Ana pipe-0.45m Ana pipe-0.70m Ana pipe-0.95m 53 Table 3.15: Calculated analytical parameters using turbulent flow over plate with uniform surface heat flux approach. Simulated Time: 10 hours. Pipe Length(m) Temperatures (K) 0.45 341 52404 2.68 245 0.7 338 78649 2.80 343 0.95 337 105800 2.81 436 Table 3.16: Calculated analytical parameters using turbulent flow over plate with uniform surface heat flux approach. Simulated Time: 15 hours. Pipe Length(m) Temperatures (K) 0.45 342 53258 2.63 246 0.7 339 79882 2.75 346 0.95 338 107444 2.77 439 Table 3.17: Calculated analytical parameters using turbulent flow over plate with uniform surface heat flux approach. Simulated Time: 20 hours Pipe Length(m) Temperatures (K) 0.45 344 54391 2.57 248 0.7 340 80965 2.70 347 0.95 340 109141 2.72 442 Table 3.18: Analytical convection coefficients using turbulent flow over plate with uniform surface heat flux approach. Time (hours) Values of convection coefficients (h) W/m 2 K @ 0.45m @ 0.70m @ 0.95m 5 353 318 297 10 358 321 301 15 361 324 303 20 365 326 306 54 Comparison between the numerical and analytical (turbulent flow over plate with uniform surface heat flux approach) convection coefficients for the same sections and times is shown in Fig. 3.17. The result shows that, at any time, the percentage differences between the numerical and analytical coefficients at the 0.95 m, 0.7 m and 0.45 m sections are 30%, 15% and 6% respectively. The results of Figs. 3.16 and 3.17 clearly demonstrate that the convective flow in the system studied is closer to a boundary layer flow, since the convection coefficients obtained numerically are of the same order of magnitude as the boundary layer convective coefficients calculated; validating the heat transfer results obtained by COMSOL Multiphysics. Figure 3.17: Analytical (turbulent flow over plate with uniform surface heat flux approach) and numerical convection coefficients at different sections as a function of time. 200 220 240 260 280 300 320 340 360 380 400 5 7 9 11 13 15 17 19 C o n ve ct io n c o e ff ic ie n t (W /m 2 K ) Time (hours) Num-0.45m Num-0.70m Num-0.95m Ana plate-0.45m Ana plate-0.70m Ana plate-0.95m 55 3.3 Phase Change Process Validation 3.3.1 Analytical Solution The phase change heat transfer encountered in a PCM storage is a transient, non-linear phenomenon, presenting a moving solid-liquid interface. Nonlinearity is the source of difficulties in moving boundary problems and, therefore, analytical solutions for phase change problems are only known for very simple physical geometries and systems, having simple boundary conditions. A well known analytical solution for a one- dimensional moving boundary problem, called Stefan’s problem, was first presented by Neumann (Jegadheeswaran et al., 2009). In order to validate the numerical results obtained for a transient solid-liquid phase change process analytically, this simplified 1D analytical model was studied. Figure 3.18 shows the schematic representation of this 1D problem: the melting of a semi-infinite slab of PCM (paraffin wax). For the present analysis, it is assumed that the entire paraffin wax is initially at its melting temperature (Tm), i.e., an infinitesimal increase in its temperature will result in phase change from solid to liquid. The initial and boundary conditions are similar to those used in the numerical simulation, and will be presented in the next subsection. Figure 3.18: Schematic of the 1D analytical phase change problem studied. Liquid Solid, Tm Melting front PCM Insulation Tw y x x=0 t=0 56 To analytically solve the presented problem, the following assumptions, defining Stefan’s problem, were made: 1. Initially, the solid PCM and the applied wall temperature are considered to be at the PCM melting temperature Tw = Tm; 2. The temperature distribution is considered to be 1D in the x-direction; 3. Conduction is the only heat transfer mechanism in the PCM and the PCM thermophysical properties are independent of temperature. The governing transient energy equation to be solved in the liquid phase is (Massoud, 2002): Where l is the liquid thermal diffusivity, defined as The boundary conditions are given as: Where X(t) is the solid-liquid melting interface position. The heat flux coming from the liquid phase at the melting interface, as well as the magnitude of the latent heat of the PCM, control the movement of the solid-liquid interface in the system through the following melting front energy balance (Massoud, 2002): The transient temperature distribution in the liquid is obtained by solving Eqs. (3.26) to (3.29) (Massoud, 2002): 57 With the similarity variable η is defined as: and the constant is obtained from Eq. (3.29) (Massoud, 2002): The Stefan number, Stel, defined as: The position of the melting front, measured from the wall at x = 0, as a function of time, is given by (Massoud, 2002): and the melting front velocity is obtained with (Massoud, 2002): The rate of heat transferred at the melting interface, or the rate of energy stored through phase change of the solid phase, is given by (Massoud, 2002): 3.3.2 Numerical Analysis Stefan’s problem is solved using a pure conduction transient modeling approach for the numerical simulation. Figure 3.19 presents a schematic of the numerical model used. The geometry is made up of a 2D slab of PCM (paraffin wax), 0.28 m in width and 0.1 m in height. The PCM is initially at it melting temperature of 313K, therefore any increase in temperature will result in melting. The specified thermal boundary conditions are such 58 that only on the left is heated (350K), while the horizontal surfaces as well as the right surface are insulated. Figure 3.19: Schematic of the 2D numerical model used to validate phase change heat transfer. In order to account for the latent heat of fusion at the melting interface, the heat capacity method is applied as presented by Eq. (2.24). The results of the numerical simulations will effectively capture the melting interface and the temperature distribution in the PCM. A mesh convergence study was performed to ensure that the right element size was used in order for the numerical results to be mesh independent. The simulations were performed for four different element sizes. The element sizes, number of elements, degree of freedom and simulation time needed by COMSOL Multiphysics are presented in Table 3.19. Table 3.19: Variation of element sizes for the PCM Element size (mm) Number of elements Number of DOF solved for Simulation time (seconds) 30 116 263 9 10 738 1553 34 3 8052 16363 118 0.6 224612 450493 784 Insulation y x x=0 t=0 59 Since the modified heat capacity method renders the system of equations to be solved non-linear, the transient simulation needs to be performed using small time steps less than 350 seconds, for a total simulated time of 75 000 seconds (21 hours), to ensure that any calculated temperature increase is never large enough to make a particular region of the system numerically jump over the phase change region, thereby passing the melting process and rendering the simulation useless. Figure 3.20 shows the linear region P-P, only 0.07 m in length, at which the numerical temperature profile will be compared to the analytical profile. The numerical simulation obtained after 21 hours of simulated time is also presented on this figure. As can be seen, after 21 simulated hours, only a quarter of the system has seen a temperature increase, justifying the use of this geometry to simulate a semi-infinite system. Figure 3.20: Simulated temperature plot for the numerical model. Simulated Time: 21 hours The convergence study performed for the phase change process in terms of the temperature distribution along section P-P for a simulated time of 3 hours is presented in Fig. 3.21. The numerical solution converges to the proper solution as the element size is reduced from 30 mm to 3 mm; further reduction of the element size to 0.6 mm no longer P P t=0 x=0 B 60 results in any noticeable change in the numerical solution. For that reason, element size of 3 mm is used in the simulation. Figure 3.21: Temperature profile for various element sizes, along section P-P. Simulated Time: 3 hours. The temperature distribution along section P-P obtained after 0.8, 3, 6, 10 and 16 hours is presented in Fig. 3.22. The result shows that 0.015, 0.028, 0.039, 0.051 and 0.063 m of PCM were fully melted during these 5 time periods respectively. This demonstrates proper simulation of heat conduction along the length of the liquid PCM. 310 315 320 325 330 335 340 345 350 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 T e m p e ra tu re ( K ) Distance from the wall (m) 30 mm 10 mm 3 mm 0.6 mm 313 314 315 316 317 318 0.025 0.03 0.035 0.04 0.045 T e m p e ra tu re ( K ) Distance from the wall(m) 61 Figure 3.22: Temperature distribution as a function of the distance from the wall, along section P-P, for various simulated times. The temperature history for the phase change process experience at point B of Fig. 3.20 for a PCM melting over a 3 K temperature range (313 to 316 K), obtained after 21 simulated hours is presented in Fig. 3.23. Melting starts at 313K and proceeds until 316K at 16 hours where an inflexion point is observed. This temperature behavior of the PCM with time is as expected for a material melting between temperatures of 313 and 316 K. 310 315 320 325 330 335 340 345 350 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 T e m p e ra tu re ( K ) Distance from the wall (m) Num-0.8hours Num-3hours Num-6hours Num-10hours Num-16hours 62 Figure 3.23: Temperature as a function of time, at point B, for a PCM with a melting temperature range of 313K-316K. Simulated Time: 21 hours. 3.3.3 Validation The comparison between the numerical and analytical temperature profiles is shown in Fig. 3.24. At 0.8 hour, the melting behavior obtained numerically (full line) corresponds exactly to the behavior predicted analytically (dash line). However, the numerical results for longer periods of time tend to under predict the total heat flux resulting in a smaller amount of melting, and lower temperature in the liquid PCM phase. This can be explained by the fact that, for this simulation, the PCM melts over a 3K temperature range; while, for the analytical model, the PCM melts at a unique constant temperature. 310 315 320 325 330 335 340 0 5 10 15 20 T e m p e ra tu re ( K ) Time (hours) 63 Figure 3.24: Numerical and analytical temperature profile in the PCM along section P-P for various time. Figure 3.25 presents the results obtained when changing the PCM melting temperature range from 3K, to 2K and down to 1K, while comparing these numerical results to the analytical solution for 21 hours simulated time. The gap between the analytical temperature profile and the numerical ones shrinks when the melting temperature range is reduce from 3 K, to 2K, and finally to 1 K. The heat capacity ) was modified accordingly following Eq. (2.26) each time the melting temperature range was changed. This result clearly demonstrates the viability of using the modified specific heat method to simulate solid-liquid phase change. It also shows the limitation of the analytical solution when it comes to the effect of the mushy region on the temperature profile and the heat flux in the system. 310 315 320 325 330 335 340 345 350 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 T e m p e ra tu re ( K ) Distance from the wall (m) Ana-0.8 hr Num-0.8 hr Ana-3 hrs Num-3 hrs Ana-6 hrs Num-6 hrs Ana-10 hrs Num-10 hrs Ana-16 hrs Num-16 hrs 64 Figure 3.25: Numerical and analytical temperature profile along section P-P for various numerical melting temperature ranges. Simulated Time: 21 hours. Further validation of the numerical solid-liquid phase change results was performed using the melting front positions, as a function of time, calculated numerically and analytically (Eq. (3.34)), as shown in Fig. 3.26. Since a mushy region is defined in the numerical model, the melting front position is chosen as the middle of this mushy region (313 to 316 K) at any given time. Due to the large discontinuity found in both models at t = 0 s, the initial position of the melting front is calculated after 0.1 hour. The results show a good agreement between the numerical and analytical melting front position. 315 320 325 330 335 340 345 350 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 T e m p e ra tu re ( K ) Distance from the wall (m) Analytical Melting range(313K-314K) Melting range(313K-315K) Melting range(313K-316K) Decrease in melting temperature range 65 Figure 3.26: Numerical and analytical melting front position along section P-P as a function of time 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0 2 4 6 8 10 M e lt in g f ro n t p o si ti o n ( m ) Time (hours) Numerical Analytical 66 CHAPTER 4: LHESS NUMERICAL ANALYSIS: RESULTS AND DISCUSSION Having shown that COMSOL Multiphysics can adequately simulate fluid flow, heat transfer and more importantly, phase change through the analytical validation of the numerical processes presented and discussed in the previous chapters, this chapter will present the results of the numerical study performed on the cylindrical axisymmetric 2D LHESS. The simulation presented in this chapter were perform for 12 hours simulated time, which is typical of the daily charging period when the storage system is use in conjunction with a solar thermal storage system. The mesh element optimum sizes and time steps used conform to those obtained in the validation models. The results present the effects of the thermal fluid’s velocity, and the number and distribution of fins on the phase change process inside the LHESS. It is also presented a calculation of the overall energy stored and the storage system efficiency of the LHESS as a function of the number of fins used. 4.1 Calculation of Energy Stored in the LHESS From the finite element simulation results, knowing the temperature at every point in the system at every time, it is possible to calculate the total amount of energy stored, as well as the contribution from sensible and latent heat, at every instant during the charging process. This can be achieve by performing a volume integration, on the PCM domain, on the temperature multiplied by Cp,pcm(T) (as define by Eq. (2.27)) divided by the total volume of the PCM: 67 By limiting the integral on region where the final temperature is found within certain ranges, the total amount of sensible and latent heat can also be computed. In COMSOL, these volumetric integral translate into the integration of the following function to account for the total energy stored, sensible and latent respectively: Equation (4.3) accounts for the sensible energy stored in the solid phase (T < 313 K) and in the liquid phase (T > 316 K), while Eq. (4.4) only accounts for the latent heat stored when the PCM melts between temperatures of 313 and 316 K. 4.2 Effects of Thermal Fluid Velocity In order to study the effects of the thermal fluid velocity on the thermal behavior and the phase change process in the LHESS, simulations of LHESS having 0, 1, 2, 3, 6, 9, 12, 15, 16, 17 and 18 fins were performed. The analysis was performed for thermal fluid velocities in the turbulent range, having Reynolds numbers of 7000, 14000, 42000, 70000 and 140,000, which translate to thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s respectively, and mass flow rates of 0.11, 0.23, 0.69, 1.14 and 2.3 kg/s respectively. A 68 mass flow rate of 0.11 kg/s would be typical of solar storage systems, but the other mass flow rates were used to study the effect of the thermal fluid velocity on the LHESS. The initial and boundary conditions were those presented in Chapter 3 (Table 3.3) and remain the same throughout this chapter. Figures 4.1 and 4.2 present the simulated temperature distribution plots for 6 and 15 fins LHESS arrangement and for thermal fluid inlet velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. Figure 4.1: Simulated temperature distribution plots for the 6 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. More heat is transferred from the thermal fluid to the PCM and fins closer to the inlet, resulting in higher temperatures of the PCM close to the pipe wall and fin surfaces and subsequently a higher melting fraction (phase change) of the PCM. Figure 4.1 shows that the heat transferred and energy stored increases with an increase of the thermal fluid velocity, which is directly related to a similar increase in Reynolds number. Velocity= 0.05m/s Re = 7000 Velocity= 0.1m/s Re = 14000 Velocity= 0.3m/s Re = 42000 Velocity= 0.5m/s Re = 70000 Velocity= 1m/s Re = 140000 69 Figure 4.2: Simulated temperature distribution plots for the 15 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Due to the higher number of fins used in the simulation results presented in Fig. 4.2, more heat is transferred from the thermal fluid to the PCM, when compared to the simulated results presented in Fig. 4.1. This also shows that more energy is being stored as the number of fins increases for the same inlet velocity. The difference in energy stored in the LHESS can also be visualized by the PCM melting interface represented by two black contour lines, one representing a temperature of 313 K, the other, 316 K. In Fig. 4.1, both contour lines can be seen, while in Fig. 4.2, more melting of the PCM is observed since only the 316 K contour line can be seen in some of the LHESS compartments. This corresponds to higher energy storage. The entire PCM melted in the case of 15 fins with a thermal fluid inlet velocity of 1 m/s, resulting in a nearly completely charged LHESS. Appendix A presents the simulated Velocity= 0.05m/s Re = 7000 Velocity= 0.1m/s Re = 14000 Velocity= 0.3m/s Re = 42000 Velocity= 0.5m/s Re = 70000 Velocity= 1m/s Re = 140000 70 temperature distribution plots for every LHESS configuration studied and for the 5 inlet velocities considered. The variation of the energy stored in the LHESS with time for a 12 hours charging period is graphically presented in Figs. 4.3 to 4.13 for thermal fluid inlet velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s and for LHESS arrangement having 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fins. The 0 fin LHESS configuration do not show a significant increase in the energy stored as the thermal fluid velocity is increased, because of the high thermal resistance present on the PCM side of the system compared to the thermal resistance found on the thermal fluid side. But as the number of fin is increase, diminishing the thermal resistance on the PCM side, the thermal resistance on the HTF side becomes relatively bigger and starts playing a bigger role on the overall heat transfer process. As a result, the effect of increasing the HTF velocity becomes more and more pronounced, resulting in more energy stored, as the amount of fins is increased. This is shown in Table 4.1 which presents the percentage increase in total energy stored for the 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin configurations as the thermal fluid velocity is increased from 0.05 m/s to 1 m/s. The percentage Energy increase is calculated based on the following equation: Table 4.1: Percentage total stored energy increase in the LHESS as thermal fluid velocity is increase from 0.05 m/s to 1 m/s for various numbers of fins. Number of fins 0 1 2 3 4 5 6 9 12 15 18 Energy increase (%) 3 20 27 31 35 39 42 55 63 62 58 The percentage increase of stored energy reaches a maximum of 63% for the 12 fin configuration. With further addition of fins, this percentage increase value drops to 62% and 58% for the 15 and 18 fin configurations respectively. This reduction is caused by the fact that for these two configurations, with an HTF velocity of 1 m/s, the entire PCM is melted. As a result, there is no more contribution of latent heat to the total energy stored, only a much smaller contribution from sensible heat storage in the liquid PCM. This can be observed on Figs. 4.12 and 4.13. 71 Figure 4.3: Total energy stored in the 0 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s Figure 4.4: Total energy stored in the 1 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 0 1 2 3 4 5 6 7 8 9 10 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 72 Figure 4.5: Total energy stored in the 2 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. Figure 4.6: Total energy stored in the 3 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. 0 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 T o ta l E n e rf g y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 73 Figure 4.7: Total energy stored in the 4 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. Figure 4.8: Total energy stored in the 5 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. 0 5 10 15 20 25 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 0 5 10 15 20 25 30 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 74 Figure 4.9: Total energy stored in the 6 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. Figure 4.10: Total energy stored in the 9 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 0 5 10 15 20 25 30 35 40 45 50 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 75 Figure 4.11: Total energy stored in the 12 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. Figure 4.12: Total energy stored in the 15 fin LHESS arrangement as a function of time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. 0 10 20 30 40 50 60 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 0 10 20 30 40 50 60 70 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 76 Figure 4.13: Total energy stored in the 18 fin LHESS arrangement as a function of the time for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1m/s. The amount of sensible energy stored in the 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin LHESS arrangements with thermal fluid inlet velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s after 12 simulated hours is presented in Fig. 4.14. The percentage sensible energy increase for the 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin LHESS configurations as the thermal fluid velocity is increases from 0.05 m/s to 1 m/s is presented in Table 4.2. Table 4.2: Percentage sensible energy increase in the LHESS as thermal fluid velocity is increase from 0.05 m/s to 1 m/s for various numbers of fins. Number of fins 0 1 2 3 4 5 6 9 12 15 18 Energy increase (%) 2 12 16 18 19 20 21 27 58 102 120 The results shows that there is a progressive increase from 2% to 27% in the sensible heat increase percentage as the number of fins is being increased from 0 to 9 fins. For these configuration, the bulk of the PCM is has not yet fully melted, and the thermal resistance offered by the PCM still dominates the heat transfer process, making the increase in HTF 0 10 20 30 40 50 60 70 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 77 velocity less significant. The sensible heat increase percentage goes to 58%, 102% and 120% for the 12, 15 and 18 fin configurations respectively. This is as a result of a large fraction of the PCM being melted, so a larger fraction of sensible heat is stored in the liquid PCM. Also, for these large number of fin configurations, the overall thermal resistance offered by the PCM/fins is smaller, making the effect of increasing the HTF velocity more noticeable. Figure 4.14: Sensible heat energy stored after 12 hours in 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin LHESS arrangements for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. Figure 4.15 presents the latent energy stored after 12 simulated hours with thermal fluid inlet velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s, for 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin LHESS configurations. For all the LHESS fin configurations presented, the latent energy increases as the thermal fluid velocity is being increased. The effect of the thermal fluid 0 5 10 15 20 25 0 1 2 3 4 5 6 9 12 15 18 S e n si b le H e a t S to re d ( M J) Number of fins 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 78 on the latent energy stored also shows that it increases with the addition of fins. This can be seen in Table 4.3, which presents the percentage latent energy increase for the 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin configurations as the thermal fluid velocity is increased from 0.05 m/s to 1 m/s. Table 4.4: Percentage latent energy increase in the LHESS as thermal fluid velocity is increase from 0.05 m/s to 1 m/s for various numbers of fins. Number of fins 0 1 2 3 4 5 6 9 12 15 18 Energy increase (%) 4 29 38 43 50 56 62 64 65 45 35 Figure 4.15: Latent heat energy stored after 12 hours in 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin LHESS arrangements for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. 0 5 10 15 20 25 30 35 40 0 1 2 3 4 5 9 12 15 18 L a te n t H e a t S to re d ( K J) Number of fins 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 79 The percentage latent energy increase for 0 to 12 fin configurations shows a progressive increase, because of overall decrease of the thermal resistance on the PCM side with an increase number of fins. Further addition of fins reduces the available PCM storage volume which reduces the amount of latent energy stored. Also, a larger amount of PCM melts at lower velocities in the 15 and 18 fin configurations, which results in a reduction in the stored latent energy increase percentage. Figure 4.16: Total energy stored after 12 hours in 0, 1, 2, 3, 4, 5, 6, 9, 12, 15 and 18 fin LHESS arrangements for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1m/s. Figures 4.16 presents the total energy stored in the 0, 1, 2, 3, 4, 5, 6, 9, 12, 15, and 18 fin LHESS arrangements for thermal fluid inlet velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s after 12 hours of charging. The total energy is the combination of the sensible and latent energy stored in the LHESS. For all LHESS fin configurations, after 12 hours of charging, the total energy stored increases as the thermal fluid velocity increases. For 0 10 20 30 40 50 60 70 0 1 2 3 4 5 6 9 12 15 18 T o ta l E n e rg y S to re d ( M J) Number of fins 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 80 thermal fluid velocity increase from 0.05 m/s to 1 m/s, the percentage increases in the total energy stored for the 0, 1, 2, 3, 4, 5, 6, 9, 12, 15, and 18 fin LHESS configurations are presented in Table 4.1. Figure 4.17 presents the PCM melting fraction in the 0 to 18 fin configurations for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s after 12 simulated hours. The temperature of 314.5 K (mid-temperature of the mushy region) was selected as the solid- liquid transition temperature for the purpose of this calculation. The results show that for any particular LHESS fin configuration, the melting fraction increases with an increase in the thermal fluid velocity: this can be seen with the 12 fin configuration in which the melting fractions for 0.05, 0.1, 0.3, 0.5 and 1 m/s thermal fluid velocities are 0.55, 0.68, 0.85, 0.9 and 1 respectively. The result also shows that complete melting (melting fraction of 1) was achieved for the 17 fin configuration for thermal fluid velocities of 1, 0.5 and 0.3 m/s. Also, the results show that there is, for each HTF velocity, a nearly linear increase in the melting fraction with an increase in the number of fins until complete melting is achieved. Figure 4.17: Melting fraction as a function of the number of fins for thermal fluid velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 2 4 6 8 10 12 14 16 18 M e lt in g f ra ct io n Number of fins 0.05m/s 0.1m/s 0.3m/s 0.5m/s 1m/s 81 4.3 Effects of addition of Fins This section describes the effect of the number of fins on the charging of the LHESS. The analysis was performed based on the results obtained from the numerical simulations carried out on 0 to 27 fin LHESS configurations. The simulated time used was 12 hours and in order to check the influence of the thermal fluid velocity on the total energy stored and subsequently the overall system efficiency, thermal fluid velocities of 0.05 m/s and 0.5 m/s were used, these velocities translate into mass flow rates of 0.11 kg/s and 1.14 kg/s respectively. From a practical point of view, considering the length and diameter of the pipe used for the LHESS (Table 2.1), it would be more practical to use a fluid velocity of 0.05 m/s (which translate to 20 seconds of thermal fluid flow through the pipe) in an actual system. 4.3.1 Effects of Addition of Fins (HTF velocity of 0.05m/s) Figure 4.18 presents the simulated temperature distribution plots for the 1, 5, 10, 14 and 18 fin LHESS configurations for thermal fluid velocity of 0.05 m/s. As the number of fins is increased, more energy is being stored in the PCM. This is shown in the temperature distribution in the PCM and the behavior of the melting front (313 K and 316 K) represented by the two black contour lines. This double melting front is fully present in the 1 and 5 fin configurations, partially present in the 10 fin configuration and as the number of fins is increased to 14, melting progresses to a temperature of 316 K and the melting front becomes single. As the number of fins is further increased to 18, a partial single melting front is observed. Appendix B presents the simulated temperature distribution plots for every LHESS configurations studied with thermal velocity of 0.05 m/s. 82 Figure 4.18: Simulated temperature distribution plots for the 1, 5, 10, 14 and 18 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.05 m/s. The variation of total energy stored with time for 0 to 27 fin LHESS configurations is presented in Fig. 4.19. The energy stored increases steadily with time as the number of fins added to the LHESS increases. However, there is a larger energy storage increase with the addition of fins when the system as fewer fins than when it has more. After the storage system was charged for 12 hours, the minimum and maximum total energy stored are 3.6 MJ and 39.8 MJ for 0 fin and 27 fins LHESS configurations respectively. One fin Five fins Ten fins Fourteen fins Eighteen fins 83 Figure 4.19: Total energy stored in some of the 0 to 27 fin LHESS arrangements as a function of time. Thermal fluid velocity of 0.05 m/s. Figure 4.20 presents the amount of sensible and latent energy stored in the 0 to 27 fin LHESS configurations for a thermal fluid velocity of 0.05 m/s. There is a progressive increase in the latent energy stored from 1.5 MJ for 0 fin to a maximum of 28 MJ for 27 fins. This is due to the increase in the melting rate of the PCM as a result of the addition of fins. The sensible energy stored increases progressively from 2 MJ for 0 fin to 10.7 MJ for 9 fins. It becomes almost constant with further addition of fins. This is because of low thermal fluid velocity resulting in a low heat transfer rate to the PCM, therefore the bulk of the energy transferred is used to melt the PCM. 0 5 10 15 20 25 30 35 40 45 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0 fin 1 fin 2 fins 3 fins 4 fins 5 fins 6 fins 7 fins 8 fins 9 fins 10 fins 11 fins 12 fins 13 fins 14 fins 15 fins 16 fins 17 fins 18 fins 20 fins 22 fins 24 fins 27 fins 84 Figure 4.20: Sensible and latent heat energy stored in the 0 to 27 fin LHESS arrangement after 12 hours of charging. Thermal fluid velocity of 0.05 m/s. 4.3.2 Effects of Addition of Fins (HTF velocity of 0.5 m/s) The simulated temperature distribution plots for the 1, 5, 10, 14 and 18 fin LHESS configurations for a thermal fluid velocity of 0.5 m/s is presented in Fig. 4.21. Again, more energy is being stored in the PCM as the number of fin is increased. This is shown in the temperature distribution in the PCM and the behavior of the melting front (313 K and 316 K) represented by the two black contour lines. This double melting front is present in the 1 and 5 fin configurations, but the melting progresses to a temperature of 316 K as the number of fin is increased to 10 and 14. As the number of fins is further increased to 18, the entire PCM is melted; the melting front disappeared totally. Appendix C presents the simulated temperature distribution plots for every LHESS configurations studied with an inlet velocity of 0.5 m/s. 0 5 10 15 20 25 30 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 20 22 24 27 E n e rg y S to re d ( M J) Number of fins Sensible Heat Latent Heat 85 Figure 4.21: Simulated temperature distribution plots for the 1, 5, 10, 14 and 18 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.5 m/s. The variation of energy stored with time in the 0 to 27 fin LHESS configurations for thermal fluid velocity of 0.5 m/s for a simulated time of 12 hours is presented in Fig. 4.22. The energy stored increases steadily with time, as the number of fins added to the LHESS is increased, the minimum and maximum total energy stored are 3.6 MJ and 39.4 MJ for 0 and 27 fins respectively. The maximum amount of stored sensible, latent and total energy in the LHESS as a function of the number of fins (available PCM volume) is presented in Table 4.4. The amount of sensible and latent energy stored in the 0 to 27 fin LHESS configurations is presented in Fig. 4.23. The latent energy stored increases steadily from 1.6 MJ for 0 fin to 35.8 MJ for 15 fins and then starts to decrease as the number of fins added to the storage system increases, until it get to 33.6 MJ for 27 fins. One fine Five fins Ten fins Fourteen fins Eighteen fins 86 Figure 4.22: Total energy stored in some of the 0 to 27 fin LHESS arrangements as a function of time. Thermal fluid velocity of 0.5 m/s. This decrease in the amount of stored latent energy is caused by the reduction in available PCM volume coming from the addition of the fins to the LHESS as shown in Table 4.4. The stored sensible energy increases steadily from 2.1 MJ for 0 fin to 12 MJ for 7 fins. Then, the increase rate reduces greatly between the 7 and 12 fin configurations creating a plateau, this translates to higher latent energy storage. Between 13 and 27 fins, the rate of sensible energy stored increases again more rapidly from 15.7 MJ to 22.9 MJ because energy is being stored in larger quantities in the liquid PCM. 0 10 20 30 40 50 60 0 2 4 6 8 10 12 T o ta l E n e rg y S to re d ( M J) Time (Hours) 0 fin 1 fin 2 fins 3 fins 4 fins 5 fins 6 fins 7 fins 8 fins 9 fins 10 fins 11 fins 12 fins 13 fins 14 fins 15 fins 16 fins 17 fins 18 fins 20 fins 22 fins 24 fins 27 fins 87 Table 4.4: Maximum stored sensible, latent and total energy in the LHESS as a function of the number of fins. Number of fins PCM Volume (m 3 ) Sensible Energy (MJ) Latent Energy (MJ) Total Energy (MJ) 0 0.286 2.12 1.61 3.73 1 0.284 4.19 4.31 8.50 2 0.283 6.21 6.91 13.12 3 0.281 8.19 9.49 17.68 4 0.280 9.93 12.03 21.96 5 0.278 11.23 14.85 26.08 6 0.277 11.97 17.99 29.96 7 0.276 12.34 21.27 33.61 8 0.274 12.64 24.53 37.17 9 0.273 12.92 27.58 40.50 10 0.271 13.21 30.30 43.51 11 0.270 13.72 32.63 46.35 12 0.268 14.62 34.25 48.87 13 0.267 15.78 35.20 50.98 14 0.266 17.09 35.63 52.72 15 0.264 18.18 35.81 53.99 16 0.263 19.35 35.75 55.10 17 0.261 20.31 35.57 55.88 18 0.260 20.99 35.38 56.37 20 0.257 22.01 34.98 56.99 22 0.254 22.50 34.60 57.10 24 0.251 22.84 34.21 57.05 27 0.247 22.94 33.63 56.57 Figure 4.24 present the results of the energy stored in the LHESS for 0 to 27 fin configurations and for thermal fluid velocities of 0.05 m/s and 0.5 m/s. This is the combination of sensible and latent energy stored for each thermal fluid velocity. The result for the thermal fluid velocity of 0.05 m/s shows progressive increase in energy as the numbers of fins are increase, from 3.6 MJ for 0 fin to 39.7 MJ for 27 fins. The energy stored for a thermal fluid velocity of 0.5 m/s also shows a progressive increase from 3.7 MJ for 0 fin to 57 MJ for 24 fins, but the energy stored decrease to 56.5 MJ for 27 fins because the volume of the PCM have been extremely reduce due to the addition of fins to the LHESS. 88 Figure 4.23: Sensible and Latent heat energy stored in the 0 to 27 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.5 m/s. Figure 4.24: Total energy stored in some of the 0 to 27 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.5 m/s. 0 5 10 15 20 25 30 35 40 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 20 22 24 27 E n e rg y S to re d ( M J) Number of fins Sensible Heat Latent Heat 0 10 20 30 40 50 60 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 20 22 24 27 T o ta l E n e rg y S to re d ( M J) Number of fins 0.05m/s 0.5m/s 89 4.4 Energy Stored versus Maximum Storage Capacity and Energy Storage System Efficiency Figure 4.25 present the total energy stored after 12 hours of simulated time by using thermal fluid velocities of 0.05 m/s and 0.5 m/s compared with the maximum energy storage capacity of the LHESS, for the 0 to 27 fin configurations. The maximum energy storage capacity is calculated using Eq. (4.6) and is based on the assumption that the entire PCM is now at the HTF inlet temperature of 350 K: The highest maximum energy storage capacity is 66.8 MJ for the 0 fin configuration, while the lowest maximum energy storage capacity obtained with the 27 fin LHESS configuration is 57.8 MJ. The linear decrease in the maximum energy storage capacity is due to the reduction in the available storage volume of the PCM as the number of fins added to storage system takes part of the available PCM volume. The maximum energy stored in the LHESS for thermal fluid velocities of 0.05 m/s and 0.5 m/s after 12 hours of charging are 39.8 MJ and 57 MJ respectively, and were obtained with 27 and 24 fins configurations. The energy storage rate and the amount of PCM melted at any given time in the 0.05 m/s is smaller than that of 0.5 m/s thermal fluid velocity, therefore it took the 27 fin configuration to achieve the maximum energy storage with an HTF velocity of 0.05 m/s. The 0.5 m/s thermal fluid velocity has a higher storage rate and all the PCM is melted with 15 fins. This explains the slow increase in total stored energy for every configuration having more than 15 fins. The reduction in available PCM volume also explains the slow increase in the total stored energy of the fully charged LHESS configurations for an HTF velocity of 0.5 m/s. 90 Figure 4.25: Comparison of total energy stored with maximum storage capacity for the 0 to 27 fin LHESS arrangement after 12 hours of charging and with thermal fluid velocities of 0.05 and 0.5 m/s. Equation (4.7) was use to calculate the LHESS efficiency. This equation depends on the number of fins used, the charging time and the thermal fluid velocity: Figure 4.26 present the LHESS efficiencies for the 0 to 27 fin configurations for thermal fluid velocities of 0.05 m/s and 0.5 m/s after 12 hours of charging. The efficiencies increase with the addition of fins. The 0 fin configuration has the lowest efficiencies of 5.4 and 5.5% for thermal fluid velocities of 0.05 and 0.5 m/s respectively, while the highest efficiencies of 68.9% and 97.9% occur with 27 fins for thermal fluid velocities of 0.05 and 0.5 m/s respectively. 0 10 20 30 40 50 60 70 0 5 10 15 20 25 T o ta l E n e rg y S to re d ( M J) Number of fins Maximum storage capacity Energy stored(0.5m/s) Energy stored(0.05m/s) 91 Figure 4.25 shows that for LHESS used with a thermal fluid velocity of 0.05 m/s, the total energy stored with 15 fins is 35.5 MJ, which corresponds to an efficiency of 57.5%. Increasing the number of fins by 44%, that is from 15 to 27, only increases the total energy from 35.5 MJ to 39.8 MJ (16.6% increase). In the case of the LHESS running with a thermal fluid velocity of 0.5 m/s, Figs. 4.25 and 4.26 shows that the total energy stored and the efficiency increases almost linearly from 3.7 MJ for 0 fin to 54 MJ for 15 fins, which corresponds to an efficiency of 87.4%. Increasing the number of fins in the system by 44%, that is from 15 to 27, increases the total energy stored by only 2.5 MJ (4.4% increase). From this discussion, using a LHESS configuration with 15 fins is a good trade-off between the total amounts of energy stored, efficiency for multiple HTF velocity and the material and manufacturing expenses needed to continuously add more fins to this system. Figure 4.26: LHESS storage efficiencies for 0 to 27 fins LHESS arrangements after 12 hours of charging with thermal fluid velocities of 0.05 and 0.5 m/s. 0 10 20 30 40 50 60 70 80 90 100 0 5 10 15 20 25 L H E S S E ff ic ie n cy (% ) Number of fins Efficiency(0.5m/s) Efficiency(0.05m/s) 92 CHAPTER 5: CONCLUSION Fluid flow, heat transfer and phase change processes were studied numerically by using the finite element method in order to simulate the charging process encountered in a cylindrical LHESS. These numerical studies were validated analytically. The effect of the addition of fins to the LHESS and the effect of the thermal fluid velocity on the storage process were studied. In this research, COMSOL Multiphysics has been successfully used to implement and model the phase change heat transfer process encountered in the cylindrical storage device studied. In order to account for the large amount of energy stored in the LHESS as a result of the large value of the latent heat of fusion of the PCM (174 kJ/kg), the melting process was simulated by modifying the PCM specific heat capacity over the melting temperature range of 3K used, resulting in an effective specific heat of 60.5 kJ/kg∙K. 5.1 Conclusion The effect of the thermal fluid velocity on the LHESS was studied. The studies presented results of various thermal fluid velocities as they affect different fin configurations of the LHESS. The effect of the thermal fluid velocity on the LHESS was studied for thermal fluid velocity ranging from 0.05 m/s to 1 m/s, and for 0 to 18 fin configurations. Results shows that after 12 hours of charging, there is a gradual increase of sensible energy stored as the thermal fluid velocity is being increased. This is the combination of sensible energy stored in the solid and liquid phase of the PCM. It increases from a minimum of 2 MJ for 0 fin and thermal fluid velocity of 0.05 m/s to a maximum of 23 MJ for 18 fins and a thermal fluid velocity of 1 m/s. 93 Due to the large amount of energy stored during melting of the PCM, the stored latent energy increases rapidly as the thermal fluid velocity and the number of fins are increased. The latent energy increases from 1.5 MJ for 0 fin and 0.05 m/s thermal fluid velocity to a maximum of 36 MJ for 12 fins and 1 m/s thermal fluid velocity. However, the amount of latent heat stored reduces to 35 MJ as the number of fins is increased to 18. This reduction in the latent energy is as a result of the drastic reduction in the available PCM volume due to the addition of fins. In general, the total energy stored which is the combination of sensible and latent energy increases steadily with an increase in thermal fluid velocity and number of fins. The results of the effect of addition of fins, from 0 to 27, for the LHESS with thermal fluid velocities of 0.05m/s and 0.5 m/s, shows that the total energy stored in the LHESS after 12 hours increases with the addition of fins. For thermal fluid velocity of 0.05 m/s, the total energy stored increases progressively from 3.6 MJ for 0 fin to 39.8 MJ for 27 fins, which correspond to a 91% energy increase due to the addition of fins. The amount of stored energy for a thermal fluid velocity of 0.5 m/s increases from 3.7 MJ for 0 fin to 57 MJ for 24 fins. The amount of energy stored decreases to 56.5 MJ for the 27 fin configuration because the available volume of PCM has been extremely reduced due to the addition of fins to the LHESS. The effect of the thermal fluid velocity and addition of fins on the LHESS shows that for configurations with a lower number of fins, increasing the thermal fluid velocity do not result in an appreciable increase in the total energy stored. This is due to the excessively large thermal resistance on the PCM/fin side compared to the convection thermal resistance on the HTF side. As the number of fins was increased, the thermal resistance on the PCM/fins side is reduce, resulting in the convection thermal resistance on the HTF side playing a relatively larger role in the overall heat transfer process. In that case, changing the HTF velocity has a bigger effect in the total amount of energy stored in the LHESS. In choosing the right LHESS fin configuration, the amount of energy stored as a function of the number of fins and the fluid velocity must be looked at, always keeping in mind the potential material and manufacturing costs stemming from the addition of fins. For 94 LHESS used with a thermal fluid velocity of 0.05 m/s, the total energy stored with 15 fins is 35.5 MJ, which corresponds to an efficiency of 57.5%. Increasing the number of fins by 44%, that is from 15 to 27, only increases the total energy from 35.5 MJ to 39.8 MJ (16.6% increase). In the case of the LHESS running with a thermal fluid velocity of 0.5 m/s, the total energy stored and the efficiency increases almost linearly from 3.7 MJ for 0 fin to 54 MJ for 15 fins, which corresponds to an efficiency of 87.4%. Increasing the number of fins in the system by 44%, that is from 15 to 27, increases the total energy stored by only 2.5 MJ (4.4% increase). From this discussion, using a LHESS configuration with 15 fins is a good trade-off between the total amounts of energy stored, efficiency for multiple HTF velocity and the material and manufacturing expenses needed to continuously add more fins to this system. Analysis of the LHESS system efficiencies with respect to the maximum storage capacity shows that the efficiencies for the 0.05 m/s and 0.5 m/s thermal fluid velocities for 27 fins configuration after 12 hours of energy stored are 68.9% and 97.9% respectively, but for the recommended number of fins of 15, for thermal fluid velocities of 0.05 m/s and 0.5 m/s, the system efficiencies are 57.5% and 87.4% respectively. This result shows that increasing the thermal fluid velocity increases the efficiency of the LHESS. However, based on actual solar domestic hot water system and the physical dimensions of the storage device, a thermal fluid velocity of 0.05 m/s would be more representative of the actual velocities encountered in such system. 5.2 Recommendations Further research on this thesis can be extended in the following areas: 1. In this research, the type and configuration of fins used are the internal radial fins that were equally spaced along the cylindrical storage device. In future studies, more complex fin configurations and arrangements can be use to achieve a better result in terms of the storage rate. 2. In order to further validate the numerical results, experimental validation, based on the present geometry, should be carried out. 95 3. This study is limited only to the charging process of the LHESS. Further work can be carried out on the discharging process of the energy stored in the LHESS, since this process might behave in an unexpected way and require a different optimal configuration. 4. Lastly, based on the geometry and the working temperature range of the present LHESS design, it is mainly suitable for solar energy application, but effort can be made to design a LHESS that will have a broader range of working applications: some of the examples include passive storage in building, thermal protection of food and electronic devices and thermal comfort in vehicles and spacecraft 96 REFERENCES ABHAT A. (1983) “Low temperature latent heat thermal energy storage: heat storage materials”, Solar Energy Journal, Vol. 30, pp. 313–332. AGYENIM F., PHILIP E. and MERVYN S. (2009) “A comparison of heat transfer enhancement in a medium temperature thermal energy storage heat exchanger using fins”, Solar Energy, Vol. 83, pp. 1509–1520. AGYENIM F., HEWITT N., EAMES P and SMYTH M. (2008) “Numerical and experimental development of medium temperature thermal energy storage (Erythritol) system for the hot side of LiBr/H2O air conditioning applications”, Proceedings of the World Renewable Energy Congress 2008. ANDREW M., MOHAMMED F., SELMAN J.R. and SAID A. (2006) “Thermal conductivity enhancement of phase change materials using a graphite matrix”, Applied Thermal Engineering, Vol. 26, pp. 1652–1661. ANICA T. (2005) “Experimental and numerical investigations of heat transfer during technical grade paraffin melting and solidification in a shell and tube latent thermal energy storage unit”, Solar energy journal, Vol. 79, pp. 648-660. ASAKO Y. and FAGHRI M. (1999) “Effect of density change on melting of unfixed rectangular phase-change material under low gravity environment”, Numerical Heat Transfer Part A, Vol. 36, pp. 825. ASSIS E., KATSMAN L. and LETAN R. (2007) “Numerical and experimental of melting in a spherical shell”, International Journal of Heat and Mass Transfer, Vol. 50, pp. 1790-1804. ATUL S., TYAGI V. V., CHEN C.R and BUDDHI D. (2009) “Review on thermal energy storage with phase change materials and applications”, Renewable and Sustainable Energy Reviews, Vol. 13, pp. 318–345. BENJAMIN J.J., DAWEI S., SHANKUR K. and SURESH V. G. (2006) “Experimental and Numerical study of melting in a cylinder”, International Journal of Heat and Mass Transfer, Vol. 49, pp. 2724-2738. BERTRAND O. and BINET B. (1999) “Melting driven by natural convection. A comparison exercise: first results”, International Journal of Thermal Sciences, Vol. 38, pp. 5-26. BRENT, A.D. VOLLER V.R. and REID K.J. (1988) “Enthalpy-porosity technique for modeling convection–diffusion phase change: application to the melting of a pure metal”, Numerical Heat Transfer, Vol. 13, pp. 297–318. 97 CHOI J.C. and KIM S.D. (1992) “Heat transfer characteristics of a latent heat storage system using MgCl2_6H2O”, Energy, Vol. 17, pp. 1153–64. CRANK J. and GUPTA R. (1995) “Isothermal migration in two dimensions”, International Journal of Heat and Mass Transfer, Vol. 18, pp. 1101–17. COMSOL Multiphysics User’s Guide, (2008) Version 3.5a. DEBORAH A. K. and MICHAEL K. J. (2005) “Introduction to Thermal and Fluid Engineering”, John Wiley and Sons, New Jersey, pp. 420–436. DEMIREL Y. and PAKSOY H.O. (1993) “Thermal analysis of heat storage materials”, Thermo Chemical Acta, Vol. 213, pp. 211-221. DINCER, I. (2002) “Thermal energy storage systems as a key technology in energy conservation”, International Journal of Energy Research, Vol. 26, pp. 567-588. DUTT S.S., HIROAKI K. and KAZUNOBU S. (2004) “Phase Change Materials for Low Temperature Solar Thermal Applications”, Research Report of Faculty of Engineering, MIE University, Vol. 29, pp. 31-64. EAMES I. W., KAMEL T. and ADRE F. (2002) “Freezing and melting of water in spherical enclosures of the type used in thermal (ice) storage systems”, Applied Thermal Engineering, Vol. 22, pp. 733–745. ESAM M.A. (2004) “Phase change process with free convection in a circular enclosure: Numerical simulations”, Computer and Fluids Journal, Vol. 33, pp.1335-1348. ERMIS K., EREK A, and DINCER I. (2007) “Heat transfer analysis of phase change process in a finned-tube thermal energy storage system using artificial neural network”, International Journal of Heat and Mass Transfer, Vol. 50, pp.3163–75. EVA G., HARALD M. and STEFAN H. (2007) “Modeling of subcooling and solidification of phase change materials”, Modeling and Simulation in Materials Science and Engineering, Vol. 15, pp. 879–892. FUKUSAKO S. and YAMADA M. (1999) “Melting heat transfer inside ducts and over external bodies”, Experimental Thermal and Fluid Science, Vol. 19, pp. 93-117. GHASEMI B. and MOLKI M. (1999) “Melting of unfixed solid square cavities”, International Journal of Heat and Fluid Flow, Vol. 20, pp. 446-452. GO Z., LIU H, and LI Y. (2004) “Thermal energy recovery of air-conditioning system heat recovery system calculation and phase change material development”, Applied Thermal Engineering Journal, Vol. 24, pp. 2511-2526. 98 GONG Z. and MUJUMDAR A. (1998) “A finite element model for convection- dominated melting and solidification problem”, International Journal of Numerical Method of Heat Fluid Flow, Vol. 4, pp. 393. GROULX, D. and LACROIX, M. (2006) “Study of close contact melting of ice from a sliding heated flat plate”, International Journal of Heat and Mass Transfer, Vol. 49, pp. 4407-4416. GROULX, D. and LACROIX, M. (2007) “Study of the effect of convection on close contact melting of high Prandtl number substances”, International Journal of Heat and Mass Transfer, Vol. 46, pp. 213-220. HALE, D.V., HOOVER, M.J., and O'NEILL, M.J., (1991), “Phase change materials Hand Book", Marshal Space Flight Center, Alabama, pp. 43-46. HAMADA Y., OHTSU W. and FUKAI J. (2003) “Thermal response in thermal energy storage material around heat transfer tubes: effect of additives on heat transfer rates”, Solar Energy, Vol. 51, pp. 317–28. HASNAIN, S. (1998) “Review on sustainable thermal energy storage technologies, Part I: heat storage materials and techniques”, Energy Conservation and Management, Vol. 39, pp.1127-1138. HENDRA R, HAMDANI T.M.I., and MAHLIA H.H. M. (2005) “Thermal and melting heat transfer characteristics in a latent heat storage system using mikro”, Applied Thermal Engineering, Vol. 25, pp. 1503–1515. HORBANIUC B, DUMITRASCUA G and POPESCUB A. (1999) “Mathematical models for the study of solidification within a longitudinally finned heat pipe latent heat thermal storage system”, Energy Conversion & Management, Vol. 40, pp.1765–74. INCROPERA F.P.,DEWITT D.P.,BERGMAN T.L. and LAVINE A.S. (2007) “Fundamental of Heat and Mass transfer”, John Wiley and Sons, 6th Edition, New Jersey, pp. 96–200. JEGADHEESWARAN S. and SANJAY D. P (2009) “Performance enhancement in latent heat thermal storage system: A review", Renewable and Sustainable Energy Reviews, Vol. 13, pp. 2225–2244. KHILLARKA D.B., GONG Z. X. and MUJUMDAR A.S.(2000) ”Melting of a phase change material in concentric horizontal annuli arbitrary cross section”, Applied Thermal Engineering Journal, Vol. 20 , pp. 893-912. KIM C. J. and KAVIANY M. (1992) “A numerical method for phase-change problems with convection and diffusion”, International Journal of Heat Mass Transfer, Vol. 35, pp. 457–467. 99 LAMBERG P. (2003) “Mathematical modeling and experimental investigation of melting and solidification in a finned phase change material storage”, Doctoral Dissertation, Department of Mechanical Engineering, Helsinki University of Technology, Espoo, Finland. LAMBERG P. and SIREN K. (2003) “Analytical model for melting in a semi-infinite PCM storage with an internal fin”, Heat Mass Transfer, Vol. 39, pp. 167–76. LAMBERG P., LEHTINIEMI R. and HENELL A.M. (2004) “Numerical and experimental investigation of melting and freezing processes in phase change material storage”, International Journal of Thermal Sciences, Vol. 43, pp. 277–87. LANE, G. (1983) "Latent heat materials", Volume 1, CRC Press, Boca Raton, Florida. LANE, G.A. (1999) "Hand book of thermal design, Chapter 1: Phase change thermal storage materials", Edited by C. Guyer, McGraw Hill Book Co., pp. 67-74. LAUARDINI V.J. (1981) “Heat transfer in cold climates”, Van Nostrand, New York, pp. 98-103. LAZARIDAS, S. (1990) “A numerical solution of the multidimensional solidification (or melting) problem”, International Journal of Heat Mass Transfer, Vol. 13, pp.1459–77. LIU, Z. and CHUNG, D.D.L. (2005) "Calorimetric evaluation of phase change materials use as thermal interface materials", Thermo Chemical Acta, Vol. 366, pp. 135-147. MASSOUD K. (2002) “Principles of Heat Transfer”, John Wiley and Sons, New York, pp. 306–311. MBAYE M, and BILGEN E. (2001) “Phase change process by natural convection– diffusion in rectangular enclosures”, Heat Mass Transfer, Vol. 37, pp. 35-44. MOORE F.E. and BAYAZITOGLU Y. (1984) “Melting within a spherical enclosure”, Journal of Heat Transfer, Vol. 104, pp. 19–23. NATERER G.F. (2003) “Heat Transfer in Single and Multiphase Systems”, CRC Press LLC, Boca Raton, Florida, pp. 356-389 NAYAK K. C., SAHA, S. K. and DUTTA P. (2006) “A numerical model for heat sinks with PCMs and thermal conductivity enhancers”, International Journal of Heat and Mass Transfer, Vol. 49, pp. 1833-1844. REGIN A. F., SOLANKI S. C. and SAINI J.S. (2006) “Latent heat thermal energy storage using cylindrical capsule: Numerical and experimental investigations”, Renewable Energy, Vol. 31, pp. 2025–2041. 100 ROLPH III W. and BATH E.K.J. (1982) “An efficient algorithm for analysis of nonlinear heat transfer with phase change”, International Journal of Numerical Methods Eng., Vol. 18, pp.119–34. ROSEN M.A. and DINCER I. (2003) “Energy methods for assessing and comparing thermal storage systems”, International Journal for Energy Research, Vol. 27, pp. 415- 430. RO S. T., LEE J. S. and SUH J. S. (1990) “Experimental study of the melting process in a spherical enclosure”, Proceedings of ICHMT 90, Yugoslavia, 1990, pp. 54-61. SAITOH T. and MOON J. H. (1998) “An experimental study for combined close-contact and natural convection melting in a spherical capsule”, International Journal of Fluid Mechanics Research, Vol. 25, pp. 220-229. SASAGUCHI K. and TAKEO H. (1994) “Effect of the orientation of a finned surface on the melting of frozen porous media”, Heat and Mass Transfer, Vol. 37, pp. 13–26. SIMPSON J.E. and GARIMELLA, S. V. (1998) “An investigation of solution, thermal and flow fields in unidirectional alloy solidification”, International Journal of Heat Mass Transfer, Vol. 41, pp. 2485–2502. SUWONDO A., MANSOORI G. and HIRAN S. (1994) “Characterization of alkanes and paraffin waxes for application as phase change energy storage medium”, Energy Sources, Vol. 16, pp. 117-128. SZIMMAT J. (2002) “Numerical simulation of solidification process in enclosures” Heat Mass Transfer, Vol. 38, pp. 72-79. URO S. (2004) “An experimental study of enhanced heat transfer in rectangular PCM thermal storage”, International Journal of Heat and Mass Transfer, Vol. 47 pp. 2841– 2847. VELRAJ R., SEENIRAJ R.V., HAFNER B., FABER C. and SCHWARZER K. (1997) “Experimental analysis and numerical modeling of inward solidification on a finned vertical tube for a latent heat storage unit”, Solar Energy, Vol. 60, pp. 281–290. VELRAJ R, SEENIRAJ RV, HAFNER B, FABER C, and SWHARZER K. (1999) “Heat transfer enhancement in a latent heat storage system”, Solar Energy, Vol. 65, pp.171–80. VISWANATH R. and JALURIA Y. (1993) “A comparison of different solution methodologies for melting and solidification problems in enclosures”, Numerical Heat Transfer Part B: Fundamental, Vol. 24, pp. 77–105. WANG Y. AMIRI A. and VAFAI K. (1999) “An experimental investigation of the melting process in a rectangular enclosure”, International Journal of Heat and Mass Transfer, Vol. 42, pp. 3639-3642. 101 WOLFF F. and VISKANTA R. (1988) “Solidification of a pure metal at a vertical wall in the presence of liquid superheat”, International Journal of Heat Mass Transfer, Vol. 31, pp. 1735–1744. YIMER B., and ADAMI M. (1997) “Parametric study of phase change material energy storage systems for space applications”, Energy Conversion and Management Journal, Vol. 38, pp. 253-262. YINGQIU Z., YINPING Z., YI J, and YANBING K. (1999) “Thermal storage and heat transfer in phase change material outside a circular tube with axial variation of the heat transfer fluid temperature”, International Journal of Solar Energy Engineering, Vol. 121, pp. 145–149. ZALBA B., JOSE M. and LUISA F. (2003) “Review on thermal energy storage with phase change: materials, heat transfer analysis and applications”, Applied Thermal Engineering, Vol. 23, pp. 251–283. ZHANG H. PRASAD V. and MOALLEMI M. K. (1996) “Numerical algorithm using multizone adaptive grid generation for multiphase transport processes with moving and free boundaries”, Numerical Heat Transfer Part B: Fundamental, Vol. 29, pp. 399–421. ZHANG Y. and FAGHRI A. (1996) “Heat transfer enhancement in latent heat thermal energy storage system by using the internally finned tube”, International Journal of Heat Mass Transfer, Vol. 39, pp. 3165–73. ZIVKOVIC, B. and FUJII, I. (2001) “Analyses of isothermal phase change of phase change material within rectangular and cylindrical containers”, Solar Energy, Vol. 70, pp. 51–61. 102 APPENDIX A: Temperature Plots as a Function of the Thermal Fluid Velocity This appendix shows the temperature plots as a function of the thermal fluid velocity. The results display are for thermal fluid inlet velocities of 0.05, 0.1, 0.3, 0.5 and 1 m/s as it affects the LHESS fin configurations of 0, 1, 2, 3, 4, 5, 6, 9, 12, 15, 16, 17 and 18 fins. The 2 black contour lines are representing the limit of the melting mushy region from 313K to 316K. Figure A.1: Simulated temperature distribution plots for the 0 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 103 Figure A.2: Simulated temperature distribution plots for the 1 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Figure A.3: Simulated temperature distribution plots for the 2 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 104 Figure A.4: Simulated temperature distribution plots for the 3 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Figure A.5: Simulated temperature distribution plots for the 4 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 105 Figure A.6: Simulated temperature distribution plots for the 5 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Figure A.7: Simulated temperature distribution plots for the 6 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Velocity= 0.1m/s Re = 14000 Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.05m/s Re = 7000 Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 106 Figure A.8: Simulated temperature distribution plots for the 9 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Figure A.9: Simulated temperature distribution plots for the 12 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 107 Figure A.10: Simulated temperature distribution plots for the 15 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Figure A.11: Simulated temperature distribution plots for the 16 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 108 Figure A.12: Simulated temperature distribution plots for the 17 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Figure A.13: Simulated temperature distribution plots for the 18 fin LHESS arrangement for various thermal fluid velocities after 12 hours of charging. Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 Velocity= 1m/s Re = 140000 Velocity= 0.5m/s Re = 70000 Velocity= 0.3m/s Re = 42000 Velocity= 0.1m/s Re = 14000 Velocity= 0.05m/s Re = 7000 109 APPENDIX B: Temperature Plots as a Function of Number of Fins (HTF Velocity of 0.05 m/s) This appendix shows the temperature plots as a function of the number of fins. The results display are for thermal fluid inlet velocities of 0.05 m/s as it affects the LHESS fin configurations of 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 20, 22, 24 and 27 fins. The 2 black contour lines are representing the limit of the melting mushy region from 313K to 316K. Figure B.1: Simulated temperature distribution plots for 0, 1, 2, 3 and 4 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.05 m/s. 0 fin 1 fin 3 fins 2 fins 4 fins 110 Figure B.2: Simulated temperature distribution plots for 5, 6, 7, 8 and 9 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.05 m/s. Figure B.3: Simulated temperature distribution plots for 10, 11, 12, 13 and 14 fin LHESS arrangements after 12 hours of charging. Thermal fluid vel. of 0.05 m/s. 5 fins 6 fins 7 fins 8 fins 9 fins 10 fins 11 fins 12 fins 13 fins 14 fins 111 Figure B.4: Simulated temperature distribution plots for 15, 16, 17 and 18 fin LHESS arrangements after 12 hours of charging. Thermal fluid vel. of 0.05 m/s. Figure B.5: Simulated temperature distribution plots for 20, 22, 24 and 27 fin LHESS arrangements after 12 hours of charging. Thermal fluid vel. of 0.05 m/s. 15 fins 16 fins 17 fins 18 fins 20 fins 22 fins 24 fins 27 fins 112 APPENDIX C: Temperature Plots as a Function of Number of Fins (HTF Velocity of 0.5 m/s) This appendix shows the temperature plots as a function of the number of fins. The results display are for thermal fluid inlet velocities of 0.5 m/s as it affects the LHESS fin configurations of 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 20, 22, 24 and 27 fins. The 2 black contour lines are representing the limit of the melting mushy region from 313K to 316K. Figure C.1: Simulated temperature distribution plots for 0, 1, 2, 3 and 4 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.5 m/s. 0 fin 1 fin 2 fins 3 fins 4 fins 113 Figure C.2: Simulated temperature distribution plots for 5, 6, 7, 8 and 9 fin LHESS arrangements after 12 hours of charging. Thermal fluid velocity of 0.5 m/s. Figure C.3: Simulated temperature distribution plots for 10, 11, 12, 13 and 14 fin LHESS arrangements after 12 hours of charging. Thermal fluid vel. of 0.5 m/s. 7 fins 8 fins 5 fins 6 fins 9 fins 10 fins 12 fins 11 fins 13 fins 14 fins 114 Figure C.4: Simulated temperature distribution plots for 15, 16, 17 and 18 fin LHESS arrangements after 12 hours of charging. Thermal fluid vel. of 0.5 m/s. Figure C.5: Simulated temperature distribution plots for 20, 22, 24 and 27 fin LHESS arrangements after 12 hours of charging. Thermal fluid vel. of 0.5 m/s. 15 fins 18 fins 17 fins 16 fins 20 fins 22 fins 27 fins 24 fins