Numerical Investigation of the Solidification of Nanoparticle-Based Colloidal Suspensions by Yousef M. F. El Hasadi A dissertation submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Doctor of Philosophy Auburn, Alabama August 3, 2013 Keywords: Colloids, Dendrites, Nanoparticles, Phase Change Materials, Solidification, Suspensions Copyright 2013 by Yousef M. F. El Hasadi Approved by Jay M. Khodadadi, Chair, Alumni Professor of Mechanical Engineering W. Robert Ashurst, Associate Professor of Chemical Engineering Daniel W. Mackowski, Associate Professor of Mechanical Engineering Amon J. Meir, Professor of Mathematics and Statistics ii Abstract In this dissertation, analytical and numerical approaches were used for analysis of the solidification process of nanoparticle-based colloidal suspensions, with a special emphasis to those used as nanostructured-enhanced phase change materials (NePCM). A one-dimensional freezing model based on an extended version of the Rubinstein Problem will be used to analyze the solidification process of cyclohexane-copper suspensions. The chosen diameters for the nanoparticles are 7, 5 and 2 nm. The value of the initial volume fraction of the nanoparticles was varied. The rejection rate of the particles was controlled through the value of the segregation coefficient value (1 corresponding to no rejection of particles, 0.1, 0.01, and 0.001). With no particle rejection, the expedited movement of the solid-liquid interface with respect to the pure cyclohexane as the volume fraction of the particles increases is not always guaranteed for the same cold side surface temperature. However, with particle rejection, for most cases tested the solid-liquid interface is decelerated with respect to pure cyclohexane as the volume fraction of the nanoparticles is increased, and this deceleration is more pronounced as the particle size decreases. This deceleration is attributed to solidification when the rejection of the particles is switched from thermal- to solutal-controlled solidification and due to the development of a constitutionally supercooled liquid on the liquid side of the interface. The maximum attained value of the concentration at the solid-liquid interface is decreasing as the initial concentration of the particles is increased; however, the value of the interface temperature is decreased as the ii concentration of the particles is increased. The transition segregation coefficient that is the non-dimensional parameter that controls the transition from thermal- to solutal-controlled solidification is increased with the increase of the particle's volume fraction and with the decrease of the particle size. A two-dimensional model, which includes the effect of the fluid flow, will be implemented to simulate the freezing of water-copper suspension frozen from the bottom side of a square cavity. The diameter of the particles is 5 and 2 nm, and the particles' mass fraction is 10%. The model is based on the combination of a one-fluid-mixture approach with the single- domain enthalpy porosity model for phase change and assuming a linear dependence of the liquidus and solidus temperatures of the mushy zone on the local concentration of the nanoparticles subject to a constant value of the segregation coefficient. Thermal-solutal convection and the Brownian and thermophoretic effects are taken into account. The solid-liquid interface for the colloidal suspension with 5 nm particle size was almost planar throughout the solidification process. However, for the suspension with particle size of 2 nm, the solid-liquid interface evolved from a stable planar shape to an unstable dendritic structure. This transition was attributed to the constitutional supercooling effect, whereby the rejected particles that are pushed away from the interface into the liquid zone form regions of high concentration thus leading to a lower solidus temperature. Using the same two-dimensional model, the suspension of water-copper will be solidified unidirectionally from the left vertical side investigating the effect of different parameters on the thermal-solutal convection formed during the freezing process. Initially, the flow in the melt consisted of two vortices rotating in opposite directions. However, at later times only one counter clockwise rotating cell survived. Changing the material of the particle to alumina results in a iii crystallized phase with a higher concentration of particles if it is compared to that of the solid phase resulting from freezing the copper-water colloidal suspension. Decreasing the segregation coefficient destabilize the solid-liquid interface and increases the intensity of the convection cell with respect to that of no particle rejection. At slow freezing rates, the resulting crystal phase consisted of lower particle content if it is compared to that resulted from higher freezing rate. iv Acknowledgments I wish to thank my advisor Dr. Jay M. Khodadadi for his supervision and support of my doctoral research. In addition, I wish to thank Dr. W. Robert Ashurst, Dr. Daniel W. Mackowski and Dr. Amon J. Meir of the Departments of Chemical Engineering, Mechanical Engineering, and Mathematics and Statistics, respectively, for serving as my committee members. I thank all of my friends for their unlimited support during this time. I dedicate my thesis to my mother Dimitra Damianidou which I lost her very recently. She was the person that taught me to be curious and to love science. I want to thank my father Mohamed Fakhri, My two sisters Mona, and Mariam EL Hasadi for their great support all these years. Finally, I want to thank the Rautenstrauch family for the support and the enlightenment that brought into my life. Permission to reuse copyrighted materials in compilation of this dissertation was granted by the American Society of Mechanical Engineers (ASME) is acknowledged. These materials include: 1) Y. M. F. EL Hasadi, J. M. Khodadadi, 2013, ?Numerical Simulation of the Effect of the Size of Suspensions on the Solidification Process of Nanoparticle Enhanced Phase Change Materials?, Journal of Heat Transfer, Vol. 135 (5). 2) Y. M. F. El Hasadi, J. M. Khodadadi, 2012, ? Numerical Simulation of the Effect of the Size of Nanoparticles on the Solidification Process of Nanoparticle-Enhanced v Phase Change Materials?, Proceedings of the 2012 ASME Summer Heat Transfer Conference (CD ROM), Paper HT 2012-58423, 8-12 July, Rio Grand, Puerto Rico. 3) Y. M. F. El Hasadi, J. M. Khodadadi, 2013,? Numerical Simulation of Solidification of Colloidal Suspensions Inside Differentially- Heated Cavity? , Proceedings of the ASME 2013 Summer Heat Transfer Confrence, Paper number HT2013-17594, 14-19 July, Minneapolis, Minnesota, USA. This dissertation is based upon work partially supported by the US Department of Energy under Award Number DE-SC0002470. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. vi Table of Contents Abstract ........................................................................................................................................... ii Acknowledgments.......................................................................................................................... iv List of Figures .............................................................................................................................. viii List of Tables ............................................................................................................................... xiv Chapter 1 Introduction .................................................................................................................... 1 1.1 Colloidal Solidification ......................................................................................................... 1 1.2 Theoretical Modeling of Colloidal Solidification ................................................................. 2 1.3 Nanostructure-Enhanced Phase Change Materials ............................................................... 3 1.4 Research Objectives and Methodology ................................................................................. 4 1.5 Outline of the Dissertation .................................................................................................... 5 Chapter 2 ......................................................................................................................................... 8 2.1 Basic Physical Concepts of Solidification ............................................................................ 8 2.1.1 Nucleation ....................................................................................................................... 8 2.1.2 Equilibrium temperature depression ............................................................................... 2 2.1.3 Interfacial melting........................................................................................................... 3 2.3 Experimental Investigations for Solidification of non-Dilute Colloidal Suspensions .......... 9 2.4 Nanostructured-Enhanced Phase Change Materials ........................................................... 16 2.5 Theoretical Modeling of non-Dilute Colloidal Suspensions ............................................... 19 2.6 Summary ............................................................................................................................. 21 Chapter 3 ....................................................................................................................................... 28 3.1 Mathematical Model ........................................................................................................... 28 3.1.1 Similarity solution ........................................................................................................ 33 3.2 Results and Discussion ........................................................................................................ 36 3.2.1 Development of the solid-liquid interface .................................................................... 37 3.2.2 Temperature distributions in both phases ..................................................................... 42 vii 3.2.3 Concentration distributions in both phases ................................................................... 44 3.2.4 Effects of particle properties....................................................................................... 46 3.2.5 Effects of the surface temperature ................................................................................ 47 3.2.6 Transition segregation coefficient ................................................................................ 48 3.4. Summary ............................................................................................................................ 50 Chapter 4 ....................................................................................................................................... 82 4.1 Mathematical Modeling ...................................................................................................... 82 4.1.1 Governing equations ..................................................................................................... 82 4.1.2 Initial/Boundary conditions .......................................................................................... 85 4.1.3 Mixture relations and effective thermophysical/transport properties ........................... 86 4.2 Solution Procedure .............................................................................................................. 90 4.3 Results and Discussion ........................................................................................................ 91 4.4 Summary ........................................................................................................................... 100 Chapter 5 ..................................................................................................................................... 120 5.1 Introduction: ...................................................................................................................... 120 5.2 Results and Discussion ...................................................................................................... 122 5.2.1 Development of the liquid fraction and concentration profiles: ................................. 122 5.2.2 The effect of changing the segregation coefficient (k0) .............................................. 126 5.2.3 The effect of changing the solidification rate: ............................................................ 128 5.2.4 The effect of the particle size on the thermo-solutal convection ................................ 128 5.2.4 The effect of the mass fraction of the particles: ......................................................... 130 5.3 Summary: .......................................................................................................................... 132 Chapter 6 ..................................................................................................................................... 150 6.1 Concluding Remarks: ........................................................................................................ 150 6.2 Suggestions for Future Work: ........................................................................................... 151 Bibliography ............................................................................................................................... 153 viii List of Figures Figure 1.1 Variety of the solid-liquid interface morphologies for kaolinite and montmorillonite colloidal suspensions solidified from the bottom for different identified mass fractions (Peppin et al., 2007)???. 6 Figure 1.2 (a) Particle rejection form a planar interface (b) particle rejection form dendritic interface, and (c) particle engulfment by planar interface??????????????????????????? 7. Figure 2. 1 Variation of the critical nucleation radius with the degree of supercooling for water. ............ 22 Figure 2. 2 Single particle suspended in a liquid being approached by the rising solid-liquid interface. ... 23 Figure 2. 3Variation of the critical velocity with the radius of the particle for the water-gold suspension for G = 1000 K/m. ....................................................................................................................................... 24 Figure 2. 4 Variation of the critical velocity with the temperature gradient for the water-gold suspension for R = 10-6 m. ............................................................................................................................................. 25 Figure 2. 5 SEM micrograph of the final microstructure and evolution of the ice front morphology for water-alumina suspensions (Deville, 2007b). ............................................................................................. 26 Figure 2. 6 3D view of the final microstructure of frozen water-alumina suspension, showing the formation of the ice lens (Lasalle et al., 2012). ........................................................................................... 27 Figure 3. 1 The Le and St numbers with 0? and pd for cyclohexane-based colloids of (a, b) CuO, (c, d) Al2O3 and (e, f) Cu. ..................................................................................................................................... 55 Figure 3. 2 Variation of *T? with 0? and pd for cyclohexane-based colloids of (a) CuO, (b) Al2O3 and (c) Cu .......................................................................................................................................................... 56 Figure3.3 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 7 nm) for various values of 0? (no particle rejection at the interface) .............................................. 57 ix Figure 3.4 Variation of the solid-liquid interface location with time for the case of copper nanaoparticles ( pd = 5 nm) for various values of 0? (no particle rejection at the interface) .............................................. 58 Figure 3.5 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 2 nm) for various values of 0? (no particle rejection at the interface) .............................................. 59 Figure 3.6 Comparison among the variations of the solid-liquid interface with time for different copper nanoparticle sizes and 0? = 0.01 (no particle rejection at the interface) .................................................... 60 Figure 3.7 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 2 nm) and various values of 0? for 0T = 7.5 oC and surfT = 2.5 oC (no particle rejection at the interface) ..................................................................................................................................................... 61 Figure 3.8 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 2 nm) and various values of 0? for 0T = 7.5 oC and surfT = 0.5 oC (no particle rejection at the interface) ..................................................................................................................................................... 62 Figure 3.9 Variation of the nondimensional supercooling with the 0? for different values of surfT and pd = 2 nm ......................................................................................................................................................... 63 Figure 3.10 Variation of ? with surfT for different values of 0? and pd ................................................. 64 Figure 3.11 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 7 nm) and various values of 0? (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001 ................................... 65 Figure 3. 12 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 2 nm) and various values of 0? (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001 ???????... 66 Figure 3.13 Variation of the solid-liquid interface location with time for the case of copper nanoparticles with 0? = 0.01 and different rejection rates for (a) pd = 7 nm and (b) pd = 2 nm .................................... 67 Figure 3.14 Variation of temperature with *x for the case of copper nanoparticles ( pd = 7 nm) with *s = 0.002 and various values of 0? for (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001 ...................................... 68 x Figure 3.15 Variation of temperature with *x for the case of copper nanoparticles ( pd = 2 nm) with *s = 0.002 and various values of 0? for (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001 ...................................... 69 Figure 3.16 Variation of temperature with *x when *s = 0.002 for different ? values with 0? = 0.01 for (a) pd = 7 nm and (b) pd = 2 nm ............................................................................................................... 70 Figure 3.17 Variation of the interface temperature iT* with ? for different values of pd ...................... 71 Figure 3.18 Variation of *LT and *equT with *x for the case of pd = 7 nm when *s = 0.002 and 0? = 0.01 for (a) ? = 1.0, (b) ? = 0.1, (c) ? = 0.01 and (d) ? = 0.001. .................................................................... 72 Figure 3.19 Variation of *LT and *equT with *x for the case of pd = 2 nm when *s = 0.002, 0? = 0.01 for (a) ? = 1.0, (b) ? = 0.1, (c) ? = 0.01 and (d) ? = 0.001. .......................................................................... 73 Figure 3.20 Variation of concentration with *x for copper nanoparticles ( pd = 7 nm) with *s = 0.002 and various values of 0? for (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001. ............................................... 74 Figure 3.21 Variation of concentration with *x for copper nanoparticles ( pd = 2 nm) with *s = 0.002 and various values of 0? for (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001. ............................................... 75 Figure 3.22 Variation of concentration with *x with *s = 0.002 for different ? values for the case of 0? = 0.01: (a) pd = 7 nm and (b) pd = 2 nm ............................................................................................ 76 Figure 3. 23 Variation of the solid-liquid interface location with time for three materials for the case of 0? = 0.05 and pd = 2 nm: (a) ? = 1.0 and (b) ? = 0.001. ........................................................................ 77 Figure 3.24 Variation of the solid-liquid interface location with time for the case of 0? = 0.01, pd = 2 nm and for various *surfT values: (a) ? = 1.0, (b) ? = 0.1, (c) ? = 0.01 and (d) ? = 0.001. ............................ 78 Figure 3.25 Variation of temperature with *x for *s = 0.002 for the case of 0? = 0.01, pd = 2 nm and for various *surfT values: (a) ? = 1.0, (b) ? = 0.1, (c) ? = 0.01 and (d) ? = 0.001. ................................... 79 xi Figure 3.26 Variation of the expedited freezing parameter epsilon(?) with the segregation coefficient (?). .................................................................................................................................................................... 80 Figure 3.27 Comparison between the predicted freezing times based on the 1-D model and experimental data of Fan and Khodadadi (2011)??????????????????????????? . 81 Figure 4.1 Geometry of the physical model .............................................................................................. 104 Figure 4.2 Grid independence results showing (a) temperature, (b) liquid fraction and (c) mass fraction of the nanoparticles at x = 0.5 H for different grids (dp = 5 nm) at two time instants. ................................. 106 Figure 4.3 Comparison between the results of the current model and Hannoun et al. (2005): (a) stream function at t = 100 s Hannoun et al. (2005), (b) stream function at t = 100 s (current model), (c) stream function at t = 200 s Hannoun et al. (2005), (d) stream function for t = 200 s (current model) and (e) instantaneous positions of the liquid-solid interfaces ............................................................................. 106 Figure 4.4 Development of the liquid fraction field (? ) at different time instants: (a) 10 s, (b) 100 s, (c) 500 s and (d) 1,000 s for an initial mass concentration field of 10% and dp = 5 nm. ............................... 107 Figure 4.5 Development of the liquid fraction field (? ) at different time instants: (a) 10 s, (b) 100 s, (c) 500 s and (d) 1,000 s for an initial mass concentration field of 10% and dp = 2 nm. ............................... 108 Figure 4.6 Development of the nanoparticle concentration field ( w? ) at different time instants: (a) 10 s, (b) 100 s, (c) 500 s and (d) 1,000 s for an initial mass concentration field of 10% and dp = 5 nm. .......... 109 Figure 4.7 Development of the nanoparticle concentration field ( w? ) at different time instants: (a) 10 s, (b) 100 s, (c) 500 s and (d) 1,000 s for an initial mass concentration field of 10% and dp = 2 nm ........... 110 Figure 4.8 Isometric view of the nanoparticle concentration field ( w? ) at 1,000 s for an initial mass concentration field of 10% and dp = 2 nm. ............................................................................................... 111 Figure 4.9 Velocity vectors in the vicinity of the liquid-solid interface ( 8.0/3.0 ?? Hx ) at 1,000 s for an initial mass concentration field of 10% and dp = 2 nm. ....................................................................... 112 Figure 4.10 Development of the liquidus temperature at different time instants: (a) t = 10 s, (b) t = 100 s, (c) t = 500 s and (d) t = 1000 s for an initial mass concentration field of 10% and dp = 5 nm. ................ 113 xii Figure 4.11 Development of solidus temperature at different time instants (a) t = 10 s, (b) t = 100 s, (c) t = 500 s and (d) t = 1000 s for an initial concentration field of 10% and dp = 5 nm. ................................. 114 Figure 4.12 Development of the liquidus temperature at different time instants: (a) t = 10 s (b) t = 100 s (c) t = 500 s and (d) t = 1000 s for an initial mass concentration field of 10% and dp = 2 nm. ................ 115 Figure 4.13 Development of solidus temperature at different time instants (a) t = 10 s, (b) t = 100 s, (c) t = 500 s and (d) t = 1000 s for an initial concentration field of 10% and dp = 2 nm. ................................... 116 Figure 4.14 Transient development of the nanoparticle concentration profiles at x = 0.5 H and along the y- axis for dp = 5 nm. ..................................................................................................................................... 117 Figure 4.15 Transient development of the concentration profiles along the y-axis for dp = 2 nm and at (a) x = 0.0043 m (b) x = 0.005 m and (c) x = 0.0055 m. ................................................................................ 118 Figure 4.16 Solid-liquid fraction fields at t = 1000 s for: (a) w? = 0, (b) w? = 10%, dp = 5 nm and (c) w? = 10%, dp = 2 nm. ........................................................................................................................................ 119 Figure 5.1 Development of the flow field (shown by the streamlines) superimposed on the contours of the liquid fraction for Ra = 7.44x104, and copper nanoparticles of dp = 2 nm at (a) t = 10 s, (b) t = 100 s, (c) t = 500 s, and (d) t = 1000 s. ....................................................................................................................... 134 Figure 5.2 Development of the concentration field for with Ra = 7.44x104, and copper nanoparticles of dp = 2 nm at different time instants (a) t = 10 s, (b) t= 100 s, (c) t = 500 s, and (d) t = 1000 s. .................... 135 Figure 5.3 Superimposed contours of the velocity streamlines and the liquid fraction at t = 1000 s for (a) water-copper (Ra = 7.44x104) and dp = 2 nm and (b) water-alumina suspensions (Ra = 4.36x104) and dp = 2 nm. ...................................................................................................................................................... 136 Figure 5.4 Concentration field at t = 1000s for (a) water-copper, Ra = 7.44x104, and dp = 2 nm and (b) for water-alumina suspensions, Ra = 4.36x104. .............................................................................................. 137 Figure 2. 53Figure 5.5 Variation of the vertical velocity (uy(m/s)) profiles at y = 0.5H and along the x- axis at t = 1000 s, and dp = 2 nm. .............................................................................................................. 138 Figure 5.6 Variation of the concentration profiles at y = 0.5H and along the x-axis at t = 1000 s and dp = 2 nm for different particle types. .................................................................................................................. 139 xiii Figure 5.7 Superimposed contours of streamlines and liquid fraction, for Ra = 7.44x104 and dp = 2 nm at t = 1000 s for (a) k0 = 1.0, (b) k0 = 0.5, (c) k0 = 0.3, and (d) k0 = 0.1. ................................................... 140 Figure 5.8 Variation of the non-dimensional maximum vorticity strength with k0, at t = 1000 s for dp = 2 nm. ............................................................................................................................................................ 141 Figure 5.9 Superimposed contours of streamlines and liquid fraction at t = 1000 s for (a) CT = 268 K, Ra =7.44x104 and for dp = 2 nm, and (b) CT = 271 K, Ra = 2.91x104. ........................................................ 142 Figure 5.10 Concentration field at t = 1000 s, and dp = 2 nm (a) CT = 268 K, Ra = 7.44x104 and (b) CT = 271 K, Ra = 2.91x104. ............................................................................................................................... 143 Figure 2. 59Figure 5.11 Superimposed contours of streamlines and liquid fraction, for dp = 5 nm, Ra = 7.44x104, and CT = 268 K for (a) t= 10 s, (b) t = 100 s, (c) t= 500 s, and (d) t= 1000 s. ........................ 144 Figure 5.12 Variation of the non-dimensional maximum vorticity strength with time for different particle sizes ........................................................................................................................................................... 145 Figure 5.13 Contours of the concentration field, for dp = 5 nm, Ra = 7.44x104, and CT = 268 K for (a) t= 10 s, (b) t = 100 s, (c) t= 500 s, and (d) t= 1000 s. .................................................................................... 146 Figure 5.14 Superimposed contours of streamlines and liquid fraction at t = 1000 s for CT = 268 K, and dp = 2 nm (a) w? = 0.05, and (b) w? = 0.1, Ra =7.44x104. ........................................................................... 147 Figure 5.15 Variation of the non-dimensional maximum vorticity strength with time for different mass fractions of the particles. ........................................................................................................................... 148 Figure 5.16 Concentration field contours at t = 1000 s for CT = 268 K, and dp = 2 nm (a) w? = 0.05, and (b) w? = 0.1, Ra =7.44x104. ..................................................................................................................... 149 xiv List of Tables xv Table 3.2: Thermophysical properties of the nanoparticles ........................................................................ 53 Table 3.3: Solid-liquid interface locations for *t = 1000, 0? = 0.01 and different values of ? .................. 54 Table 4.1 Thermophysical and transport properties for the solvent and the nanoparticles ...................... 101 Table 4.2 Values of the dimensionless groupings at the onset of solidification (t = 0) ........................... 102 Table 4.3 Liquidus slopes and segregation coefficient used in the current study .................................... 103 1 Chapter 1 Introduction In this chapter, the motivation for the research conducted on phase change of colloidal suspensions a description of the goals of the project, and the methods used will be discussed. In the final section, the structure of the dissertation is outlined. 1.1 Colloidal Solidification Colloidal suspensions are multi-component liquid systems containing fine size particles ranging from 10-9 to 10-6 m, i.e. covering the 1- nanometer to 1- micron range. These multi- component materials are widely encountered in nature (blood and milk), as well in industrial applications such as ceramic suspensions, and nanofluids (colloidal suspensions that exhibit high thermal conductivity). Colloidal solidification is a process that occurs in nature for example the freezing of seawater and saturated soil, and in several industrial applications, such as freeze-casing of materials, purification of water, building nanorods, and cryopreservation of blood cells. In all these processes, controlling the solidification process is a decisive factor to achieve the desired results in terms of distributions of the particles. However, the mechanics of the colloidal solidification process are not understood. What has made colloidal solidification difficult to understand is the complex nature of the interaction between the suspended particles with the solid-liquid interface. Consequently, different morphologies can be obtained by changing slightly one of the controlling parameters. For instance, the variety of interfaces separating the 2 Distinct layers during freezing of two different colloidal suspensions (kaolinite, that is a type of clay, which is a mixture of alumina and oxide silicon particles in water and motmorillonite (that is a hydrated sodium calcium aluminum magnesium silicate hydroxide), crystallized from the bottom for different weight fractions of particles are shown in Figure 1.1 (Peppin et al., 2007). Three distinct morphologies of the solid-liquid interface were developed, dendritic (as in Fig. 1.1 (a), (d) and (e)), bands (Fig. 1.1 (b) and (c)), and finally hexagonal (Fig. 1.1 (f)). The variety of morphologies observed in Figure 1.1 suggests how rich and unpredictable the crystallization of colloidal suspension is. A sound theoretical knowledge is needed to thoroughly understand the nature of those morphologies and accurately predict their occurrence. Based on experimental observations, there are three possibilities for the evolution of the particles during their interaction with a solid-liquid interface. These particles might be rejected, engulfed, or trapped by the crystallized phase as show in Figure 1.2. 1.2 Theoretical Modeling of Colloidal Solidification Traditionally, there have been two approaches used to simulate the freezing of colloids, i.e.: - Single particle interactions with the solid-liquid interface, - Considering the colloidal suspension as a binary mixture, and applying the methods that have been used to simulate solidification of binary alloys. The single particle interaction with the solid-liquid interface was first studied by Uhlmann et al. (1964) who concluded that the particle would be rejected by the interface, if the velocity of the interface is lower than a critical velocity. However, the particles will be engulfed if the velocity of the interface is higher than the critical velocity. They found that the critical velocity is inversely proportional to the particle size. They showed that for small size particles, 3 the critical velocity depended strongly on the surface energies between the solid-liquid, particle- liquid, and particle-solid phases. Additional investigations were conducted to refine the understanding of the critical velocity and the parameters that affect its value by Kober and Rau (1985), Shangguan et al. (1989), and Rample and Worster (1999). However, colloidal suspensions consist of a large number of particles and using the one particle model would be computationally expensive. Moreover, the interaction between the particles plays a crucial role in the diffusion of the particles and such complexity was ignored in the single-particle models. The experiments of Peppin et al. (2008) suggest that colloidal suspensions consisting of fine particles exhibit similar behaviors as those of binary mixtures upon crystallization. These include creation of dendritic structures due to the constitutional supercooling effect. Upon considering these similarities, Peppin et al. (2008) developed a one-dimensional model of solidification based on the classical Stefan problem. Furthermore, by using a stability analysis for the modified Stefan problem, they showed that as the particle size decreases, the solid-liquid interface becomes destabilized as in Peppin et al. (2007). Finally, according to Elliot and Peppin (2011), as particles decrease in size, their potential to be captured by the crystallizing phase reduces significantly. 1.3 Nanostructure-Enhanced Phase Change Materials Nanostructure-enhanced phase change materials (NePCM) introduced by Khodadadi and Hosseinizadeh (2007) recently attracted the attention of the scientific community due to their enhanced thermal conductivity, which helps to shorten the times for charging and discharging of stored or liberated thermal energy. NePCM can be considered as colloidal suspensions with nano size structures (e.g. discrete particles, agglomerated particles, rods, wire, carbon nanotubes, graphene, etc.). Most current studies regarding the NePCM are concerned with measuring the 4 thermal conductivity of composites, as in Acme et al. (2010), Yavari et al. (2011), and Zheng et al. (2011). Few studies in the literature have focused on the performance of the NePCM during the freezing and melting cycles. Among those, Wang et al. (2011) reported a monotonic improvement of the solidification time as the volume of the nanoparticles was increased. However, the experiments of Sanousi et al. (2011), Cui et al. (2011), and Fan and Khodadadi (2012) revealed that as the volume fraction of the nanoparticles increased, the freezing or melting times did not increase monotonically. Theoretically, solidification of the NePCM was investigated by neglecting the mass transfer of the particles. This assumption represents an over- estimation, which goes against most experimental observations of colloidal solidification, and thermodynamic rules. 1.4 Research Objectives and Methodology To understand the better utilization of the NePCM in real-life applications, the solidification process associated with a multi-component system must be investigated thoroughly. This dissertation will be devoted to this purpose through completing the following research objectives: 1- Implementing a 1-D Stefan model that takes into account the mass transfer of the particles. To study the effects of rejecting the particles on the growth of the solid-liquid interface. 2- Developing a 2-D solidification model that takes into account the solutal and thermal convection effects and the variation in thermo-physical properties with the concentration of the nanoparticles. 5 3- Implementing the 2-D solidification model for water-copper NePCM, to study the effects of the particle size, volume fraction of the particles and the cold side temperature, on the development of the solid-liquid interface as well as on the variation of the concentration of the particles. 4- Investigating the effect of the rejection of the particles on the thermo-solutal convection. 1.5 Outline of the Dissertation The remainder of the dissertation consists of six chapters. A broad literature survey of the experimental and theoretical developments concerning colloidal solidification will be provided, and its connection to the thermal energy storage in Chapter 2. In Chapter 3, a one- dimensional freezing model based on the Stefan' problem will be developed, and the impact of rejecting highly conductive particles on the development of the solid-liquid interface will be investigated. A two-dimensional solidification model that accounts for the convection effects that result from the density change due to the thermal and solutal effects will be presented in Chapter 4. In Chapter 5, the effect of convection on the colloidal solidification will be investigated. Finally, in Chapter 6, concluding remarks will be presented. 6 Figure 1.1 Variety of the solid-liquid interface morphologies for kaolinite and montmorillonite colloidal suspensions solidified from the bottom for different identified mass fractions (Peppin et al., 2007). 7 (a) (b) (c) Liquid Phase , solid Phase Figure 1.2 (a) Particle rejection form a planar interface (b) particle rejection form dendritic interface, and (c) particle engulfment by planar interface. 8 Chapter 2 Literature Survey In this chapter, a literature review of solidification of colloidal suspensions will be provided. The investigations that are reviewed have focused on experimental and theoretical studies conducted for solidification of dilute suspensions, experimental work for solidification of non-dilute suspensions, phase change materials, and finally theoretical work concerning solidification of non-dilute suspensions. However, before proceeding to the literature review, the basic concepts of solidification will be outlined. 2.1 Basic Physical Concepts of Solidification 2.1.1 Nucleation During phase change, molecules of solid come out of the liquid melt through the nucleation process. Nucleation is responsible for creating the first atomic clusters of the crystalline structure that arise from random fluctuations of the liquid molecules (Dantzig and Rappaz, 2009). For the first clusters of atoms to create the first nuclei, an energy barrier ( cnhG? ) must be overcome. This free energy accounts for the free energy reduction due to the phase change, and the increase of the free energy due to the creation of the solid surface (surface effects) (Spannuth, 2010). There are two distinct nucleation processes: 1. homogeneous nucleation, and 2. heterogeneous nucleation. 1 Homogeneous nucleation occurs within a pure substance at a very slow rate. The critical radius of the nuclei represents the smallest radius (R) above which the solid cluster can grow. According to Dantzig and Rappaz (2009): slmmsc n h RTTTLRG ???? 23 4)(34 ????? (2.1) where cnhG? is the Gibbs free energy, sl? is the solid-liquid interfacial energy (J/m2), L stands for the latent heat (J/kg), s? is the density of the solid phase (kg/m3), mT is the melting temperature (K) and T is a temperature below the melting temperature. By differentiating equation (2.1) with respect to R, and by setting the result to zero, we get the value of the critical radius (Rc): )(2 TTL TR ms mslc ?? ? ? (2.2) As shown in Figure 2.1, the critical radius Rc decreases as the value of supercooling )( TTm? increases, suggesting that nucleation is easier at low temperatures, and it is the main reason that we encounter in nature the existence of a supercooled liquid, which is a liquid with a temperature below the melting temperature. Homogenous nucleation can be observed in a highly controlled environment as in a laboratory setup as Stan et al. (2009) indicated. Most of the freezing processes in nature are controlled by the heterogeneous nucleation, which is nucleation that is initiated by the presence of foreign surfaces such as the walls of the container that hold the material and any suspended particles (Dash et al., 1995). Heterogeneous 2 nucleation will only occur if the interfacial surface energy between the solid phase and the foreign surface is less than sl? . The initiation temperature of heterogeneous nucleation depends on the amount of the surfaces that surround the liquid phase of the material, as well as the materials of those surfaces (Dash et al., 2006). After nucleation, the solid-liquid interface will grow as a planar shape if the amount of supercooling is low and will be dependent upon the diffusion of heat. However, if the degree of supercooling is high, the progress of the interface will be fast, and it will depend on the kinetics of molecular attachment to the solid phase surface. In binary alloys and solidification of colloids, the liquid ahead of the solid-liquid interface can be supercooled due to its composition. This supercooling is well known in the literature as the ?constitutional supercooling.? Constitutional supercooling is the mechanism responsible for the wide range of morphologies resulting from the solidification of colloidal suspensions, as discussed by Peppin et al. (2008). 2.1.2 Equilibrium temperature depression The liquid melt can be below its equilibrium temperature (i.e. supercooling) due to the effect of the nucleation. However, there are two mechanisms that lower the equilibrium temperature (i.e. the temperature at which the liquid phase and solid phase coexist). These mechanisms are the Gibbs-Thomson effect, and the colligative effect. A brief description of these mechanisms will be given in the following paragraphs. The Gibbs-Thomson effect arises due to the contribution of the surface to the overall Gibbs free energy of the system. The increase or decrease of the free energy depends on whether 3 the curvature is concave or convex. The reduction of the equilibrium temperature due to the curvature effects is given by the following relation: LTTT s mRslmmR ????? (2.3) The depression of the equilibrium temperature due to the Gibbs-Thomson effect is particularly influential for freezing of liquids inside extremely narrow pores, or in the case of solidification of liquid droplets, or melting of solid particles. The second mechanism responsible for the depression of the freezing temperature is the colligative effect that is a result of adding solute molecules to the solvent. The depression of the equilibrium temperature is not related to the surface effect, but due to the reduction of the free energy of the suspension because of the increase of the entropy of the suspension as a result of the adding solute (Dash et al., 1995). For example, the freezing temperature of the water can be reduced to -2 oC if a certain amount of salt is added. 2.1.3 Interfacial melting Interfacial melting is a liquid layer resulting from the melting of a solid phase caused by contact with a foreign particle, due to the repulsive intermolecular interactions between the solid and the particle (Dash et al., 1995). The thickness of the liquid film resulting from the melting (Dash et al., 2006) is given by the following relation: 3131 )()6( ??? m m s H T TTLAd ?? (2.4) 4 where HA is the Hamaker constant (J). The thickness of the layer is reduced as the temperature decreases. Interfacial melting is one of the mechanisms responsible for particle rejection and engulfment during the solidification, as shown by Rample and Worster (1999). Having provided a brief overview of the physical concepts of solidification, a review of the relevant literature devoted to solidification of colloidal suspensions follows. 2.2 Experiments and Modeling for Dilute Colloidal Suspensions Corte (1962) conducted freezing experiments using particle with different materials, shapes and sizes suspended in water. The samples were frozen in a direction opposite to gravity (i.e., up from the bottom). The suspensions were frozen in cylindrical containers. The lighter particles used were made from glass while the heavier particles were made from rutile (a mineral primarily consisting of titanium oxide, i.e. TiO2). Corte (1962) observed that the lighter particles migrated greater distances from the solid-liquid interface when compared to the distances travelled by the heavier particles. The author also reported that the smaller particles migrated far from the solid-liquid interface for all the freezing rates adopted in the investigation. However, this was not the case for the bigger particles that managed to diffuse away from interface only for the low freezing rates, whereas for the high freezing rates, the particles were engulfed by the solid-liquid interface. Corte (1962) attributed this behavior to the fact that for the particle to migrate, a layer of water must exist between the particle and the crystallized phase. For small particles, the thickness of the required water layer is small while for bigger particles, a greater water flow is required to maintain the thin water layer, which is only achieved at low freezing rates. The same investigator continued the experiments (Corte, 1963) to investigate the effect of freezing and thawing cycles on the particle distribution. He observed that the particles were 5 sorted into vertical and horizontal layers with respect to that of the freezing front. He concluded that the distribution of the particles is dependent on the rate of freezing and thawing, as well as on the orientation of the freezing side and on the size of the particles, and how well dispersed in the water they are. Uhlmann et al. (1964) were the first to propose the concept of the critical velocity that formed a theoretical basis to understand when the foreign particles in a liquid will be engulfed or rejected by the liquid-solid interface. They conducted experiments with different organic/inorganic solvents (such as water, but other solvents were not specified) and particles (graphite, magnesium oxide, silt, silicon, tin, diamond, nickel, zinc, iron oxide, and silver iodide) ranging from one to several microns in diameter. They observed that at low interface velocities, in most systems the particles were pushed (rejected) by the interface even if the particle is Brownian or non-Brownian in size. However, as the velocity of the interface was increased, they observed that the particles were trapped (engulfed) in the solid phase. They found that the critical velocity depends on the particle type, and the solvents, and varies inversely with the radius of the particle. They showed that for the small size particles, the critical velocity depended strongly on the surface energies between the solid/liquid, particle/liquid and particle/solid phases. However, as the particle size increased, they found that the effect of the viscous drag force became more marked. Cisse and Bolling (1971a) conducted freezing experiments on water suspensions with copper, silicon oxide particles, hollow carbon spheres, and tungsten with the size of the particles ranging between 3 to 60 ?m. The authors intended to investigate the effect of the size of the 6 particles, the thermal conductivity, and roughness of the particle surface on the critical velocity. The suspensions that were contained in a double-walled tubing were frozen from the bottom (i.e. opposite to the direction of gravity). The interaction of the particles and the solid-liquid interface was monitored by using a light microscope. They observed one of the following three situations: (i) steady-state pushing occurred when the particles were carried hundreds of micrometers, (ii) particles could be pushed for a few micrometers before being trapped by the solid phase, or (iii) the particles were trapped by the solid phase without being pushed. They mentioned that the particles settled into an odd-shaped cone formed from the junctions of three- grain boundary grooves as the interface was pushed much farther than those that settled on the grain-boundary surface. They observed that the particles that are near the solid-liquid interface helped the nucleation of gas bubbles at the interface, which helped to trap the particles at higher interface velocities. The authors reported that the critical velocity was increased as the particle size was reduced, and the roughness of the surface of the particle is increased. They found that, as the thermal conductivity of the particles increases, their critical velocity is reduced. The authors also observed that in some cases, the solid phase that surrounds the particle is re-melted, but they could not provide any explanation of the phenomena they observed. Korber et al. (1985) investigated the interaction of the particles with the solid-liquid interface, because it plays a significant role in the cryopreservation of blood cells since most blood cells are damaged when they are exposed to a high concentration of salt in front of the solid-liquid interface. It is highly important to define the conditions that the cells will be engulfed by the crystalline phase. The authors decided to study the effect of the thermal gradient and the presence of the solute on the critical velocity of engulfment. They selected latex spheres of 5.7 and 11.9 ?m in diameter that were suspended in water and water-NaMnO4 (sodium 7 permanganate) aqueous solutions. The solution was frozen in between two glass surfaces of 20 ?m in thickness, and the interface was observed through a cryomicroscope. The authors did not maintain the velocity of the interface constant, in order to better simulate real-life scenario. The authors reported that the critical velocity is increased as the particle size is decreased. They verified their experimental results with a theoretical model similar to those of Uhlmann et al. (1964). They found that for the case of water-NaMnO4 aqueous solution, the interface was changed from a flat interface to a dendritic morphology due to the rejection of NaMnO4 (sodium permanganate) solute. However, the authors found that the critical velocity of the particles is nearly the same for the pure water and water-NaMnO4 aqueous solution for the same thermal gradient. The critical velocity was found to increase as the temperature gradient increased. Shangguan et al. (1989) investigated the interaction of the solid-liquid interface with the particle by using a steady-state model in which a balance of the repulsive and attractive forces was employed to calculate the critical velocity. Their findings were similar to that of Uhlmann et al. (1964). Rample and Worster (1999) developed a mathematical expression for the critical velocity by balancing the intermolecular and viscous forces that act on the particle. The premelted region around a single particle is shown in Figure 2.2 as represented by the premelted thickness d. They showed that a premelted liquid film will be developed around the particle and concluded that the critical velocity is inversely dependent on the size of the particle, and less sensitive to the temperature gradient applied. However, as the size of the particle decreases, the influence of the curvature of the equilibrium temperature significantly affects the shape of the liquid-solid interface. They developed the following relation for the critical velocity (Vc): 31)6( s HLA ??? ? (2.5) 8 p m sC RG T LV 41 4 1 4 9 24 3 ? ??? (2.6) Where pR represents the radius of the particle. To illustrate the effect of increasing the particle size as well as increasing the temperature gradient on the critical velocity, the data for the water-gold suspension will be used. The following values are given for the different parameters of Eq. (2.6), HA = 1.573x10-21 J, ?L 3.35x105 J/kg, ?? 8.9x10-4 Pa.s, ?s? 920 kg/m3. Figure 2.3 shows that as the particle size decreases, the critical velocity increases for a constant temperature gradient. The critical velocity increases as the temperature gradient increase for a constant particle size as shown in Figure 2.4. As shown from the two previous figures, the critical velocity is more sensitive to the change in the particle size, than the temperature gradient. However, as the particle diameter is reduced below 10-4 m, curvature effect on the melting temperature is significant. The same authors (Rample and Worster, 2001) extended their previous work to smaller size particles <10-4 m. They showed that the critical velocity is inversely proportional to the radius of the particle. However, they interestingly found that the critical velocity is not dependent on the temperature gradient imposed because the curvature effects cause a greater shift to the melting temperature than that caused by the imposed undercooling. Recently, dynamic models have been introduced to investigate the interaction of a single particle with the solid-liquid interface. In these models, the position of the particle and velocity relative to the interface are calculated based on the imbalance between the attractive drag force 9 and repulsive interfacial force, as in Azouni and Casses (1998), Catalina et al. (2000), and Garrin and Udaykumar (2005). Garvin and Udaykumar (2003a) and Garvin and Udaykumar (2003b) presented a dynamic model to investigate the engulfment or rejection of a particle approached by a solid-liquid interface. The interface and the particle were tracked using a sharp interface method. The model takes into account the difference in the thermal properties between the particle and the solvent, which helps to understand more accurately the phenomena occurring between the solid-liquid interface and the particle. No convection effects were included in the models of Garvin and Udaykumar (2003a) and Garvin and Udaykumar (2003b), however, the authors compared their results to the available experimental and analytical data from the literature. All the previous studies dealt with a single particle in the suspension. However, colloidal suspensions consist of many particles, and the interaction of the particles with each other may play a role in the solidification process. 2.3 Experimental Investigations for Solidification of non-Dilute Colloidal Suspensions In the previous section, a brief literature review was presented for the solidification process of highly dilute colloidal suspensions. Those studies were primarily concerned with the physics of the interaction of a single particle with the solid-liquid interface. In this section, an overview of the experimental work concerning the solidification of concentrated colloidal suspensions (defined as suspensions with the volume fraction of the particles higher than 1%) will be presented. Halde (1980) studied the effect of freezing on purifying water. He argued that due to the highly ordered crystalline structure of solidified water, accommodation of any foreign atom or 10 particle will not be possible, and thus most impurities will be rejected out of the solid phase. To demonstrate his argument, water with different impurities such as graphite, NaCl, and CaCO3 was frozen inside a cylindrical tube immersed in a temperature-controlled water bath while the suspension was stirred constantly. The velocity of the interface was controlled through the rate at which the container was immersed in the water bath. It was observed that increasing the stirring rate results in more purified ice than one resulting from lower stirring rates. Surprisingly, it was observed that the ice resulting from suspensions with larger particle sizes contained smaller amounts of impurities than those resulting from suspensions prepared from smaller sized particles. This contradicts the theoretical models discussed in the previous section, which explicitly showed that the small size particles would be rejected easier than the larger ones; however, the author did not give any explanation for the observation. Finally, the purification process becomes difficult as the concentration of the impurities is increased, since it causes a significant depression in freezing temperature. Another application of solidification of colloid suspensions is in investigating the interaction of the red blood cells with the ice crystal that is essential for developing appropriate protocols for cell survival. Ishiguro and Boris (1994) conducted freezing experiments with blood cells in different solutions. They used mixtures of blood cells in physiological saline, physiological saline with glycerol, and physiological saline with glycerol plus antifreeze. All the samples were frozen unidirectionally subject to different cooling rates. They observed that for the same crystallization rate, the solid-liquid interface featured different morphologies for the three different solutions: cellular, for the solution with saline, dendritic for a solution with saline with glycerol, and needle-like morphology for saline with glycerol plus antifreeze solution. The blood cells are usually trapped between the branches of the unstable interface and then engulfed 11 by the solid phase, or pushed by the solid-liquid interface. The mechanical interactions that are responsible for reshaping the cell are found to be strong for the case of saline solution, since most for the cells are damaged after the freezing. However, those interactions are inhibited by the presence of glycerol, while adding the antifreeze alters the mechanical interaction and significantly affects the morphology of the blood cell. The authors also conducted thawing experiments to test the condition of the blood cells. They reported that for the solution that contained only saline, all the cells were damaged while for the solutions that contained glycerol and antifreeze, the blood cells survived. Chang et al. (2007) conducted experimental and numerical investigations of the interaction of a single blood cell with the solid-liquid interface. They used the level-set method to follow the interface explicitly. They only solved for the concentration field and velocity in the crystalline and liquid phases since they assumed that the temperature is not affected during the solidification process because the thermal diffusivity value is much higher than the solute diffusivity in both molten and solid phases. Their model accounted for the volume change of the cell due to the expressing of water (i.e. transport of water out of the cell through the membrane) to achieve osmotic equilibrium. The long-range forces van der Waals forces on the cell were computed based on the instantaneous pair-wise summation of the dispersion forces. The drag force was calculated using the lubrication approximation. The authors used an adaptive grid that allowed them to refine the mesh in the space between the particle and the solid-liquid interface. They found that the rejection or the engulfment of the cells depends on their size. For the case of a flat interface, they observed that the cell had a higher chance of being engulfed if it is exposed for a shorter time to the high concentration of electrolytes near the interface. Under this condition, water loss was avoided and the volume of the cell did not decrease significantly. However, for the case of a dendritic interface, they found 12 that the cell was engulfed in the inter-dendritic space where the concentration of the electrolyte is higher. Freezing of food is another area in which understanding solidification of colloid suspensions is vital since it is the governing mechanism that helps to maintain desirable properties such as flavor, texture, appurtenance and mouth feel. Mashl et al. (1996) investigated the unidirectional freezing of 2% (by weight) of starch?water mixture. The suspension was contained in an optical cell, which was placed between a cold and hot side, and the interface was monitored with a light microscope. The freezing rate of the interface was controlled by the rate at which the cells moved from the hot to the cold side and the interface velocities that were achieved were in the range of 0.5 to 10 ?m/sec. The authors observed the formation of bands of high and low concentration of starch parallel to the solid-liquid interface. They reported that the morphology of the bands is altered with the change in the interface velocity. For example, at low velocities such as 0.5 and 2.0 ?m/sec, the bands tend to be more diffuse, and can easily be spotted in the microstructure. However, as the velocity of the solid-liquid interface increases to 7.5 ?m/sec, the segregation of ice is weak and the formation of the band is much more difficult to visualize. The banding phenomenon is a consequence of the change in the velocity of the solid-liquid interface and a similar microstructure was observed in binary alloy as in Elmer et al. (1994). The authors attributed the change in the interface velocity to the interaction with the particles. The starch will accumulate in front of the interface, which will result in a significant reduction of its velocity, and a decrease in the interface temperature. As the temperature of the interface decreases, the driving force for the solidification will increase, and thus the interface velocity will increase resulting in the engulfment of starch in front of the crystallizing phase. 13 However, as the liquid-solid interface advances, its temperature increases, and the velocity of the interface is reduced again, and starch is rejected out of the ice. This cycle of increasing and decreasing of the interface velocity is repeated throughout the freezing process and explains the bands observed in the frozen samples. Rejection of the particles helps to create a variety of microstructures such as cellular, dendritic, and ice lenses. The diversity of the structures created helps to develop materials of different pore sizes. The pore size depends on the microstructure that evolves from the solidification process (Deville et al., 2009a). Watanabe et al. (2001, 2002) investigated the formation of an ice lens, which can be defined as zones from which all the particles are expelled, that was formed in a perpendicular direction to the temperature gradient. The authors used a suspension of water with glass micro particles with 2 ?m in size. The suspension is directionally solidified in a specially equipped cell, capable of accommodating different cell transverse velocities. The growth of the ice is monitored with a light microscope. They observed that the formation of the ice lens grew faster as the transverse velocity of the cell increases. However, ice grew without any lens if the velocity increased beyond 2 ?m/sec. The ice was growing at a slower rate for the water glass suspension if it is compared with that of pure water. Waschkics et al. (2009 and 2011) investigated the effect of interface velocity, particle loading and size on the interface morphology. They used a water-alumina suspension that was solidified directionally from the bottom. They observed that the interface was planar for low interface velocity, with the accumulation of the particles in front of the interface. As the velocity of the interface increased, the interface took the shape of lamellars (parallel layers of a solid phase (solvent) interleaved with layers of the other phase (segregated particles)) and the particles accumulated in the space 14 between the lamellars. The lamellar interface is attributed to the constitutional supercooling that resulted from the particle pile in front of the interface. The lamellae spacing varied inversely with the interface velocity, and the spacing increases as the particle size increases. Nakagawa et al. (2010) investigated the pore size resulting from the freeze drying of a suspension of water and carbon nanotubes. The authors used carbon nanotubes with diameter size of 13-16 nm, and length of about 1-10 ?m. Two freezing methods were used to freeze the suspension: one is by contact freezing to a heat exchanger, and the second one is by immersing the freezing cell in a cryo-bath. The authors also numerically simulated the temperature history by using a one- dimensional heat conduction model. They found that the numerical results are in a good agreement with their experimental tests. They found that the pore size is related to the cooling rate and the faster the cooling rate, the smaller the pore size that developed. Furthermore, they reported that there is a pore distribution along the freezing direction, where the pores with the small size are accumulated at the bottom, while the larger pores are concentrated far from the bottom. The microstructure of the frozen sample depended on the freezing method that was used. Deville et al. (2007) investigated experimentally the two-dimensional structures that result from the freeze-drying of water-alumina suspensions that were solidified unidirectionally. The particle sizes used were 400 nm and 100 nm. The freezing cell used is surrounded with Teflon, and its shape is a cylinder of 18 mm diameter and length of 30 mm. The sample was frozen from the bottom using a cold finger. When the cold temperature at the bottom of the cell is kept constant, Scanning Electron Microscopy (SEM) images of the surface show that microstructures are lamellae and channels, but their orientation over the cross-section is entirely random. However, if the cold temperature is changed at a constant rate, which produces a higher solidification rate than in the previous case, the microstructure consists of lamellae pore 15 architecture, and there is a homogeneity and order in both perpendicular and parallel directions to the liquid-solid interface. Furthermore, some lamellae surfaces consisted of dendrite-like structures. They found that with a faster cooling rate, finer microstructure resulted. The space between the lamellars is longer for suspensions with particle sizes of 100 nm, rather than that of 400 nm. The authors clearly describe the morphologies of the evolution of the solid-liquid interface along the freezing length. They identified the following morphologies from the bottom to the top: planar where most of the particles are engulfed, columnar, columnar to lamellar, and lamellar to dendritic. For the non planar cases, the particles are rejected to the spaces between the ice crystals as shown in Figure 2.5. However, understanding the behavior of solidification of colloid suspensions through the structures resulting from the freeze-drying process may prove to be limited. Because, for example, the binder material, which has been used to keep the particles together during the sublimation process affects the viscosity of the suspension, the diffusion of the particles and the freezing temperature, and thus the final crystal structure. In-situ solidification experiments combined with X-ray radiography and tomography allows exploring in more detail the evolution of the process of solidification of colloid suspensions. X-ray can be used to monitor the solidification process of the opaque colloidal suspensions since the X-ray method depends on the absorption coefficient of the material, and thus the area of the material with significant particle concentration will appear darker than the area with less particle concentration. Another benefit of using X-rays is the high spatial resolution that will resolve the individual particles. Tomography will help to create a 3-D structure of the frozen sample showing the crystalline structure and the particle redistribution. Deville et al. (2009a, 2009b, 2009c) investigated the solidification of water-alumina suspension by using an in-situ setup with X-ray imaging and tomography. They identified that the solidification can be characterized by 16 two zones: an initial zone (where the interface velocity changes with time) and steady ones (where the interface velocity attains a constant value). The crystalline structure consists of lamellar crystals in the freezing direction, as well in the transverse directions. Lasalle et al. (2012) investigated the particle redistribution and structure formation resulting from the freezing of alumina colloids by using the same experimental setup as Deville et al. (2009a). They tried to identify the conditions under which the ice lenses are formed as shown in Figure 2.6. Ice lenses are ice crystals that form perpendicular to the direction of the heat gradient direction. They can be considered as structural defects that reduce the mechanical strength of the resulting freeze cast material. They found that the solid-liquid interface velocity is reduced as the thickness of the accumulated particles is increased. The authors found that the nucleation of the ice lenses starts as depleted zones from particles that are created above the interface. They reported that as the particle size decreases, the interface becomes unstable. Based on experiments, Peppin et al. (2008) showed that for low values of the interface velocity, the particles will be rejected away from the growing solid phase for the case of bentonite colloids, creating a denser layer of particles neighboring the solid-liquid interface. However, as the interface velocity increased, dendritic structures and cellular interface morphology developed due to the constitutional supercooling. 2.4 Nanostructured-Enhanced Phase Change Materials Recently, solidification of nano-scale colloid suspensions was identified as a potential new development in the thermal energy storage due to the introduction of the nanostructured- enhanced phase change materials (NePCM) by Khodadadi and Hosseinizadeh (2007). In the following paragraphs, a brief description of the current research on NePCM will be given, whereas a comprehensive review is also available (Khodadadi et al., 2013). Cui et al. (2011) 17 explored the feasibility of adding highly conductive nanostructures such as carbon nanotubes and carbon nanofibers to paraffin wax and soy wax (mass fractions ranging from 0 to 10 wt%). They found that the thermal conductivity at room temperature for the solid phase of the mixture was increased for both types of solvents as the mass fraction of the additive increased. They also conducted a melting experiment and found that by increasing the nanofiber content, the thermal response of the materials during extraction of thermal energy improved. However, this was not the case for the carbon nanotubes since they found that the thermal response improves up to a certain mass fraction after which it begins deteriorating. Wang et al. (2011) investigated the thermal performance of paraffin-micron size graphite flakes composites experimentally. They measured the thermal conductivity in the solid and liquid phases of the composite by using the hot disk method and found that thermal conductivity increased as the content of graphite increased. Furthermore, they observed that the values of the melting temperature and specific heat decreased as the concentration of the additive increased. In addition, they performed a solidification test with a maximum graphite flakes content of 5% by weight but they did not give any details about how the investigation was conducted. Their results indicated that the solidification time was reduced by adding graphite flakes. Sanusi et al. (2011) investigated the solidification behavior of n-tricosane (C13H28) paraffin-carbon nanofiber composites. The diameter and length of the nanofibers were 2-100 nm and 100 ?m, respectively. The mass fraction of the composite was 10 wt%. The solidification experiments were conducted for three rectangular cells with an aspect ratio of 0.5, 1 and 2, and for two heating power values of 500 and 1000 W applied on the bottom surface while the top surface of the cell was kept at a constant temperature of 5 oC. The authors reported a reduction of the solidification time with different percentages compared with pure PCM for the two aspect 18 ratio cases of 0.5 and 2. However, the authors reported an increase in the solidification time for the rectangular geometry of aspect ratio of 1 that was attributed to the change of the aspect-ratio of the geometry. However, they did not give any explanation for how a higher conductive material can solidify slower than the less conductive material. Acme et al. (2010) measured the thermal conductivity of a graphite-KNO3/NaO3 eutectic salt, which is primarily used for high- temperature thermal energy storage applications. They used the transient hot plate method to measure thermal conductivity. They found that by adding graphite to the salt, the thermal conductivity of the composite increased substantially compared to the pure salt. For example, they reported that for graphite mass fraction of 20 wt%, the thermal conductivity of the composite is about 20 Wm-1K-1. Yavari et al. (2011) measured the thermal conductivity of a 1- octadecanol-graphene composite by using a steady-state one-dimensional heat conductivity measurement method. They found that by adding 4 wt% of graphene, the thermal conductivity of the composite increased by 140% and the latent heat of the composite decreased by 15%. Furthermore, they observed that by adding more graphene to the PCM, the crystallization temperature of the composite was lowered. Wu et al. (2009) conducted simple freezing experiments on water-Al2O3 nanofluids with mass fraction of the nanoparticles ranging from 0.05 to 0.2 wt%. The samples were placed in glass tubes and then positioned in a water bath with its temperature kept constant at -15 oC. The temperature was monitored using an infrared camera. They found that by increasing the mass fraction, the freezing time reduced. Ho et al. (2009) measured the thermal conductivity, viscosity and density of an n-octadecane-Al2O3 suspension for different mass fractions of Al2O3 nanoparticles experimentally. They found that the density decreases linearly with the increase in the mass fraction of the nanoparticles. However, for thermal conductivity and viscosity, the dependence on mass fraction was nonlinear. Shin and 19 Banerjee (2011a, 2011b) reported enhancement of the heat capacity of molten salt eutectics- silica nanofluids. They attributed the enhancement of the specific heat of the colloid to one of the following mechanisms. The first mechanism is due to the increased specific heat of the nanoparticle compared to that of the bulk material, While the second mechanism is due to the layering of liquid molecules on the surface of the nanoparticle to form a semi-liquid layer. 2.5 Theoretical Modeling of non-Dilute Colloidal Suspensions In this section, a description of the analytical and non-analytical methods concerned with the modeling of solidification of non-dilute colloidal suspensions is provided. Peppin et al. (2006) developed a continuum-based model for analyzing the solidification process of hard- sphere colloidal suspensions. They showed that a Fick?s law-type transport equation is suitable to simulate the mass transport of the colloidal particles. They obtained relations for the melting temperature and the diffusion coefficient of the suspension that were nonlinear functions of the particle volume fraction. They solved a modified Stefan problem for the one-dimensional case of freezing from the bottom, assuming constant properties in the solid and liquid phases, except that for the case of the diffusion coefficient varied with the volume fraction of the particles. They did not take into account any convection or gravitational effects. Assuming that the system was under equilibrium conditions due to the slow freezing rates and all the particles were being rejected from the solid phase, the value of the segregation coefficient was set to zero. For small Lewis numbers, which represent small particles, the Brownian diffusion is strong and the particles diffuse away from the solid-liquid interface since their propagation velocity is much higher than that of the interface. As the undercooling (difference between the temperature of the cold side and the melting temperature of the solvent) increases, the concentration gradient increases at the interface and is steep enough that the gradient of the freezing temperature is greater than that of the temperature gradient. This causes the temperature ahead of the interface 20 to be colder than that of the melting temperature. This promotes the initiation of the constitutional supercooling phenomenon and the destabilization of the interface. This phenomenon is similar to what has been observed for the solidification of binary alloys. According to Peppin et al. (2006), as the Lewis number increases which represents the case of larger particles, the Brownian diffusion is weak and the particles accumulate ahead of the solid- liquid interface creating a porous layer. Peppin et al. (2007) extended their previous study to include a linear stability analysis to investigate the parameters that affect the deformation of the solid-liquid interface. They found that increasing the volume fraction of the particles and decreasing the particle radius destabilize the interface. They concluded that the mechanism responsible for the instability was the constitutional supercooling. However, the primary drawback of their model was that it did not take into account the dependence of the thermophysical and transport properties with the volume fraction and the variation of the segregation coefficient with the velocity of the interface. Elliott and Peppin (2011) recently developed a mathematical formulation for the variation of the segregation coefficient for the case of non-equilibrium solidification of colloidal suspensions. They showed that the segregation coefficient depends on the velocity of the solid-liquid interface, critical velocity and the molecular velocity of the particle. For the critical velocity, the authors used an empirical relation obtained from the literature for the case alumina-water colloidal suspensions. The segregation coefficient value was equal to zero (full rejection) at low velocities of the interface but then increased smoothly to 1 (full engulfment) as the velocity of the interface was increased. As for the effect of the size of the particle, the transition to a nonzero segregation coefficient was delayed at higher velocities as the particle diameter was lowered. They also derived a relation for the interface temperature, as well as for the non-equilibrium phase diagram. The variation of 21 the segregation coefficient with the velocity of the interface helped to clarify the banding phenomena that occurred during the experiments of Deville et al. (2009c). In these experiments, they observed zones of high and low particle concentrations in the frozen samples that indicate that the velocity of the interface is not constant but fluctuates and that the segregation coefficient does not have a constant value but changes between zero and one. 2.6 Summary A brief review of the literature for the solidification process of colloidal suspensions was presented, as well as a short overview of the literature concerning the nanostructured phase change materials (NePCM). The literature review revealed that the freezing of colloids is a complicated process and that involves different phenomena. Most investigations are concentrated on the experimental side of the problem. Especially, the morphology of the structures evolved, particle concentration profiles, and the effect of the operation conditions. The solidification of colloidal suspensions is important in many aspects of life, such as building materials with controlled properties, and structures, food industry, self-assembly processes, and health sciences. The theoretical analysis is limited to a single particle approaching the solid-liquid interface, or to solving a one-dimensional Stefan type problem for limited cases. To elucidate the experimental observations, multi-dimensional models must be developed that ccount for fluid flow. A new venue for the colloidal freezing is in the assessment of the NePCM performance as phase change material. As shown from the brief literature review presented, there is a lack of understanding of the performance of the NePCM during the charging phase (solidification), and discharging phase (melting) due to the lack of a theory that accounts for particle rejection and diffusion through the solid-liquid interface, a process that it is existed in the real-life 22 performance of the NePCM. Implementing models similar to those of colloidal freezing will help to determine the operational conditions where the NePCM can be used. Figure 2. 1 Variation of the critical nucleation radius with the degree of supercooling for water. 23 Figure 2. 2 Single particle suspended in a liquid being approached by the rising solid-liquid interface. 24 Figure 2. 3Variation of the critical velocity with the radius of the particle for the water-gold suspension for G = 1000 K/m. 25 Figure 2. 4 Variation of the critical velocity with the temperature gradient for the water-gold suspension for R = 10-6 m. 26 Figure 2. 5 SEM micrograph of the final microstructure and evolution of the ice front morphology for water-alumina suspensions (Deville, 2007b). 27 Figure 2. 6 3D view of the final microstructure of frozen water-alumina suspension, showing the formation of the ice lens (Lasalle et al., 2012). 28 Chapter 3 One-Dimensional Stefan Problem Formulation for Solidification of Colloidal Suspensions with Emphasis on Nanostructure-enhanced Phase Change Materials (NePCM) During the solidification of colloidal suspensions, a significant number of particles will be rejected from the growing solid phase. In this chapter, the effects of rejecting nanoparticles with different rates on the development of the solid-liquid interface, temperature and concentration profiles will be explored. An extended version of the Rubinstein problem will be used in the current investigation, which admits a closed-form analytical solution. The non- dimensional interface position, temperature, and concentration will be presented for the case of cyclohexane-copper suspensions with different particle sizes. It will be shown that the solidification process switches from thermally controlled to solutal controlled as the amount of the particles rejected out of the crystalline phase is increased. 3.1 Mathematical Model A thorough analysis of solidification of NePCM colloidal suspensions should account for mass transfer of the particles. Peppin et al. (2006) showed that during solidification, colloidal suspensions behave similar to a binary alloy when the particle size is a small. Another aspect of the process is the constitutional supercooling that develops due to the rejection of the nanoparticles ahead of the solid-liquid interface. The only case for which the Stefan problem 29 exhibits an analytical solution for the case of binary alloy solidification is that for the case of a semi-infinite domain cooled from one side. This is well-known as the Rubinstein problem (Crank, 1987). Voller (2008) extended the Rubinstein problem to include the an undercooled melt (i.e. the liquid phase is initially below its melting temperature). In the present work, a semi-infinite domain ( 0?x ) filled with a NePCM that is cooled from one side is considered. Initially at (t = 0), the homogeneous NePCM is at a temperature 0T and concentration 0? . The basic model and its solution method were adopted from Voller (2008), however, modifications of the thermophysical/transport properties for the NePCM are new extensions. The non-dimensional parameters used in the current study are as follows: ? ? ? ? ? ? ?? ? ? ? ? ? ?? ?? ? nL nS Le q u n npL L nL nL S L npL npS p nL nS nL npL n Lm StT L C mSt D Le m m C C C D D D l t t l s s l x x C L mTT T ? ? ?? ? ? ? ? ? ? ? ? =),(1= =,=,=,=,= =,=,=,=,= ** 0 ** 2 *** 0 0* (3.1) The superscript * represents any quantity that has no dimensions, whereas superscript + is utilized only for the dimensionless concentration. Subscripts n, L and S identify colloidal mixture properties, liquid and solid phases, respectively. Here T is temperature, mT is the melting temperature of the pure PCM or solvent, l is an appropriate length scale, L is the latent heat, pC is the specific heat, ? is the volume fraction, 0? is the initial volume concentration, s is the location of the solid-liquid interface, t is time, ? is the thermal diffusivity, D is the mass 30 diffusivity of the nanoparticles, ? is the segregation coefficient and *equT is the non-dimensional liquidus temperature. Quantities mL and mS represent liquidus and solidus slopes, respectively. The temperature at 0=x (i.e. surfT ) was lowered below the initial liquidus temperature, which is sufficient to initiate the solidification process and for the solid-liquid interface to move away from the cold surface. The governing equations in the non-dimensional form for the present ?two-region? formulation follow: Energy equation in the solid region: )(0,= ***2* *2*** tsxxTtT ?????? ? (3.2) Energy equation in the liquid melt region: )(,= ***2* *2** tsxxTtT ????? (3.3) Nanoparticle diffusion in the solid region: )(0,= ***2*2** tsxxLeDt SS ?????? ?? ?? (3.4) Nanoparticle diffusion in the liquid melt region: )(,1= ***2*2* tsxxLet LL ????? ?? ?? (3.5) The initial non-dimensional conditions of the NePCM are the following: 31 1=,= *0* ?LTT ? (3.6) The boundary conditions imposed on the boundaries of the domain are: 0=0,= **** ??? ? xatTTx s u r fS? (3.7a) ????? **0*,1 xasTTL? (3.7b) Since the Stefan problem is classified as a moving boundary problem, boundary conditions must be imposed at the solid-liquid interface. The temperature of the interface is: )(1=* iLi StT ??? ? (3.8) Where superscript i refers to the interface value and iL?? represents the nanoparticle concentration on the liquid side of the interface. Neglecting the density change across the interface, the heat balance across the solid-liquid interface is: * * * * * *** =|| dtdsxTxTC SLp ?????? (3.9) Nanoparticles are rejected from the crystalline structure due to the thermodynamics constraints and a difference of nanoparticles concentration was created between the solid and liquid sides of the interface. The concentration of the nanoparticles on the solid side of the interface is related to that of the liquid side by the following relation: iLiS ?? ??? = (3.10) Considering mass balance at the solid-liquid interface, the following relation must hold: 32 * * ** )(1=1 dtdsxLexLeD iLLLSS ? ?? ? ????? ???? (3.11) Utilizing subscript p to denote particle properties, density, specific heat, thermal conductivity, thermal diffusivity, heat of fusion, viscosity, and mass diffusivity (Brownian) diffusivity for the NePCM suspensions are calculated from the following relations: pLSLnS ????? 0,0, )(1= ?? (3.12) ppLSPLnSp CCC )())((1=)( 0,0, ????? ?? (3.13) ))(2( ))(22(= ,0, ,0,, pLSLSp pLSLSpLnS kkkk kkkkk ??? ??? ?? (3.14) LnSp LnSLnS Ck , ,, )(= ?? (3.15) Ln LL )()(1=)( 0 ??? ? (3.16) 2.50 )(1= ??? ? LnL (3.17) nLp BnL d kTD ??3= 0 (3.18) where kB, and dp are the Boltzmann constant and the diameter of the suspended particles, respectively. Utilization of a comma subscript means that either solid or liquid properties are being computed. Note that the thermophysical and transport properties are based on the initial value of the concentration field. 33 The liquidus slope is calculated from the following relation: Lv Tkm Lp mBL ? 2= (3.19) The solutions of the temperature and concentration fields along with the instantaneous position of the interface depend on ratio of mass diffusivities ( *D ), ratio of specific heats ( *pC ), segregation coefficient (? ), the Lewis number (Le), the Stefan number (St ), ratio of thermal diffusivities ( *? ), non-dimensional surface temperature ( *surfT ) and non-dimensional initial temperature ( *0T ). Upon specifying the properties of the solvent for two phases, the surface and initial temperatures, the properties of the nanoparticles, their diameter (dp), the initial volume fraction ( ?0? ) and the segregation coefficient, this problem can be solved. 3.1.1 Similarity solution The location of the moving interface can be calculated from the following relation: ** 2= ts ? (3.20) Equations (3.2) to (3.5) in combination with the initial and boundary conditions (3.6) to (3.11) admit a set of analytical solutions that are given in the following set of relations: )(0, )(e rf ) 2 (e rf )(= *** * ** * **** tsxt x TTTT s u r fis u r f ???? ? ? ? (3.21) )(0,= *** tsxiLiS ???? ??? (3.22) 34 )(,)(e r f c )2(e r f c )(= **** * *0**0* tsxt x TTTT i ??? ? (3.23) )(,)(e r f c )2(e r f c 1)(1= **** * tsxLet Lex i LL ??? ?? ??? (3.24) The values of the concentration on the liquid side of the interface iL?? , temperature of the interface iT* and ? are calculated subject to satisfying the condition of the linear variation of the liquidus concentration equation (3.8) and the modified interface balance equations (3.9) and (3.11). By substituting relation (3.20) into the interface conditions (3.9) and (3.11), one obtains: 0=)(e rf c )( )(e rf )( 2*0* * 2 ** ** * ? ? ? ? ? ???? ? ? ???? eTTeTTC is u r fi p (3.25) 0=1)()(e r f c)(1 2 ??? ?? iLLeiL LeeLe ?????? ? (3.26) Equations (3.8), (3.25) and (3.26) are solved numerically using the Newton?s method root-finding function provided by the Mathematica software (2007). The above model has the following advantages: - It is robust, easy-to-use and computationally not expensive, - It accounts for the property variations of the solid and liquid phases, - The model can incorporate the constitutional supercooling effect that is a fundamental mechanism during the solidification of colloidal suspensions. Simultaneously, the adopted model has the following disadvantages: 35 - It is one-dimensional and takes into account only conduction of heat and mass transfer of nanoparticles, however it is incapable of recovering convection effects, - It does not account for the mushy zone where the liquid and solid phases can coexist, - Temperature and concentration values at the solid-liquid interface are constant during the solidification process, - The variation of the thermophysical properties with the nanoparticle content is evaluated at the initial instant before the solidification process begins. This can be justified by the low volume fractions that have been used in this study that correspond to the dilute limit, - The suspended particles are assumed to be mono size and agglomeration of the particles subsequent to initiation of freezing is ignored. The solidification process of NePCM is characterized by two mechanisms, namely thermally-driven solidification and solute-driven solidification. In the limit of thermally-driven crystallization, the advancement of the solid phase is governed by diffusion of heat. However, for solute-driven crystallization, the growth of the solid phase is dependent on particle diffusion away from the interface. It is necessary to establish a criterion that can be used to identify which of the two mechanisms will prevail. Chiareli et al. (1994) determined that for binary alloys, a segregation coefficient defined as: )( 0msurf Lt TT m??? ?? (3.27) can be used as a measure of the importance of each mechanism. In this paper, we will name t? as the transitional segregation coefficient. According to Chiareli et al. (1994), for ? < t? diffusion of the solute will govern the process of solidification. Otherwise, thermal diffusion 36 will control the evolution of the solid phase. To realize thermally-driven solidification, Eq. (3.27) suggests that bigger particle size and greater supercooling is needed. The transitional segregation coefficient is an influential parameter that affects the performance of NePCM as shown later. 3.2 Results and Discussion Variations of the non-dimensional location of the solid-liquid interface, concentration, and temperature will be presented for cyclohexane-based NePCM. Al2O3, CuO and Cu nanoparticles with diameters of 7, 5, 3 and 2 nm that are suspended in cyclohexane (Tm equal to 6.5 oC) are considered. The thermophysical properties for the solvent and the nanoparticles are given in Tables 3.1 and 3.2, respectively. The values of the initial volume fractions considered were equal to or less than 0.01, however, higher loadings on the order of 0.05 were also studied. Variations of the Le and ?St? numbers with the initial concentration of nanoparticles and nanoparticle size are shown in Figure 3.1. These quantities are only related to the thermophysical properties. The Le number increases linearly with the initial concentration due to the dependence of the thermal diffusivity on the thermal conductivity that rises with concentration. However, the Le numbers decrease significantly as the nanoparticle size is lowered since the mass diffusivity depends inversely on the particle diameter. The slope of the ?St? number?s dependence on the initial concentration is small for the pd = 7 and 5 nm cases. However, as the particle size decreases, the ?St? number increases sharply against the initial concentration due to the dependence of the liquidus slope on the particle size. Throughout the simulations, the initial temperature of the colloid was maintained at 7.5 oC that is always above the melting temperature of the colloid, whereas the surface temperature was varied from 0.5 to 6 oC. Variations of the difference between the initial and surface 37 temperatures ( ||= **0* su rfTTT ?? ) with the initial volume fraction of the nanoparticles, as well as their size, and with surfT =4.5 oC are shown in Figure 3.2. For higher nanoparticle size, *T? is nearly constant as the value of the initial volume fraction increases. However, for smaller nanoparticle size, the values of *T? is reduced as the initial volume fraction increases. 3.2.1 Development of the solid-liquid interface The dependence of the non-dimensional location of the solid-liquid interface on the non- dimensional time for different cases will be presented here. First, for no nanoparticle rejection (suspension remains homogeneous for 0?t ) will be considered that corresponds to ? = 1 for the cyclohexane-copper suspension with pd = 7 nm, 5 nm and 2 nm, an initial temperature of 7.5 oC and cold side temperature of 4.5 oC. Note that since no particle is being rejected at the interface, according to equation (3.26), the non-dimensional concentration ( ?? ) will be equal to unity everywhere. Moreover, according to equation (3.8), the interface temperature ( iT* ) is zero. The movement of the interface is consistently accelerated as the initial nanoparticle loading is increased as shown in Figures 3.3 and 3.4 for the cases of pd = 7 nm and 5 nm, respectively. Due to the tightness of the curves, upward- and downward-pointing arrows are utilized in the legends of these figures to identify speedup and slowdown, respectively, compared to the previous loading value. These observations of expedited freezing are similar to the findings of Fan and Khodadadi (2008) who only considered a thermal conduction without including the particle size effect on the melting temperature. However, for the case of pd =2 nm, the trend associated with the monotonic speedup of the interface with particle loading is completely reversed (Figure 3.5). For this small diameter of particles, as the initial nanoparticle 38 volume fraction increases, there is a clear slowdown of the liquid-solid interface when compared with the lower values of the initial volume fraction. The movement of the interface is controlled by equation (3.20) that depends on the non-dimensional parameter ? and non-dimensional time. Parameter ? for the case of no particle rejection (? = 1) depends only on the solution of equation (3.25), that in turn varies with the values of *? , *pC and *T? . However, neither *? nor *pC depend on the size of the nanoparticles. Only *T? (difference of the dimensionless initial and surface temperatures) is a size-dependent quantity due to the dependence of the liquidus slope on the size of the nanoparticles as shown in equation (3.19) and plotted in Figure 3.2. Quantity *T? was observed to decrease as the size of the nanoparticle decreased due to the rise of the liquidus slope. Reducing the intensity of *T? that is viewed as the driving force for solidification is responsible for the slowdown that is observed in Figure 3.5. Comparing the development of the solid-liquid interface for three copper nanoparticle sizes with 0? = 0.01 (Figure 3.6), one can conclude that the solid-liquid interface for the pd = 7 nm and 5 nm nanoparticles moves nearly with the same velocity, whereas for the case of pd = 2 nm, the interface moves with a lower velocity. To investigate the effect of *T? on the movement of the solid-liquid interface for the case without rejection of nanoparticles (? =1), cyclohexane-copper suspensions with pd = 2 nm and surfT = 2.5 oC and 0.5 oC were investigated. It is observed that for the case of surfT = 2.5 oC (Figure 3.7), the interface is accelerated with an increase of the initial volume fraction if it is compared to that of pure cyclohexene. However, for the case of 0? = 0.05, the development of the solid- liquid interface is significantly slower when compared to that of all the pervious volume 39 fractions and the pure cyclohexane ( 0? =0) as well. For the case of surfT = 0.5 oC (Figure 3.8), the interface is accelerated for all the initial volume fractions of the nanoparticles studied. The variations of *T? for the cases corresponding to Figures 3.5, 3.7 and 3.8 are shown in Figure 3.9. Among these curves, that for surfT = 0.5 oC had the highest *T? value for the spectrum of the nanoparticles loadings used if it is compared to the cases of surfT = 4.5 oC and 2.5 oC. Examination of Figure 3.9 suggests that to accelerate freezing by adding highly conductive nanoparticles, a suitable non-dimensional *T? must be selected. With no particle rejection at the solid-liquid interface, the particle diameter and temperature of the cold surface determine the movement of the interface when compared to a particle-free case. Having discussed cases with no particle rejection, we now focus on cases for which mass transfer is present. The interplay of thermal- versus solutal-driven solidification is governed by the transition segregation coefficient. Figure 3.10 shows the variation of t? with the surface temperature for the cyclohexane-copper NePCM suspension for different particle sizes and initial concentrations. The value of the t? is reduced sharply as the value of the surface temperature is lowered until it attains a slower-decaying slope. The values of the t? increase as the initial concentration of the nanoparticles is raised or as the particle size is reduced. As discussed above, according to Chiareli et al. (1994) thermally-controlled solidification will dominate when the segregation coefficient is greater than the transition segregation coefficient. Thus, according to Figure 3.10, to operate an NePCM-based system for its enhanced thermal conductivity attribute (thus expedited solidification), the transition segregation coefficient should be kept small. In effect, the NePCM should be operated at low surface temperatures and/or low concentration of 40 particles. The transient development of the solid-liquid interface for the case of copper nanoparticles ( pd = 7 nm) corresponding to three values of the segregation coefficient and for different values of the initial volume fraction of the particles are shown in Figure 3.11a-c. For the case of ? = 0.1 (Fig. 3.11a), increasing the initial volume fraction of the particles results in the acceleration of the liquid-solid interface when compared to the pure cyclohexane. However, as the value of ? decreases to a value of 0.01 and 0.001 and more particles are rejected to the liquid side, the motion of the interface is decelerated as the initial volume fraction of the nanoparticles increases (Figures 3.11b and 3.11c). The interface movement for the case of 0? = 0.01 is slower than that of pure cyclohexane even though the thermal diffusivity of the former is higher than the later. The maximum value of t? for the cases considered in Figure 3.11 is 0.0046. For the case of Figure 3.11a, solidification is driven by thermal diffusion since ? > t? which explains the expedited behavior observed by increasing the initial concentration of particles and thus the thermal diffusivity, Fan and Khodadadi (2012). However, for the case of Figure 3.11b, the expedited behavior of the solid-liquid interface is not observed even though the ? > t? condition still prevails (solidification should be governed by thermal diffusion). This occurs because the liquidus and solidus temperatures are reduced with the addition of the particles, and the effect of the constitutional supercooling is increased as it will be shown in the next section. In Figure 3.11c, with ? < t? , solidification is driven by the slow diffusion of the particles. The observed declaration can be attributed to two factors. First, solidification is driven and controlled by mass diffusion of the nanoparticles that is orders of magnitude slower than the 41 diffusion of heat for the present cases of 1??Le . Second, this is due to the extended zone affected by constitutional supercooling ahead of the moving interface, as it will be shown in the next section. For the case of pd = 2 nm copper particles (Figure 3.12), the interface movement is more sensitive to the value of ? when compared to that of pd = 7 nm (Figure 3.11), where by increasing the initial volume fraction of the nanoparticles, the movement of the interface is significantly slower if it is compared to that of 0? = 0. The maximum value of t? for the cases of Figure 3.12 is 0.2029, and thus in all cases the solidification will be driven by the solutal diffusion. The significant slowdown of the interface can be attributed to the extended zone of the constitutional supercooling in front of the interface, plus that the solidification is driven entirely by the particle's diffusion. We had identified two mechanisms that are responsible for the slowdown of the interface and those are solute-driven solidification and constitutional supercooling. However, we could not identify which mechanism had the most significant effect because we could not separate each one from the other in the cases considered in this study. Increasing the thermal diffusivities in both the solid and liquid phases does not always guarantee a faster crystallization rate for smaller particle sizes as predicted by single-component solidification models. To explain this, it is argued that by reducing the particle size, the solidification will be more controlled by the solute diffusion. Experimental, observations of the interface slowdown with the addition of nanoparticles have been reported by [Sanusi et al. (2011), Fan and Khodadadi (2012), Lasalle et al. (2012), Watanabe (2002)]. Comparison of the cases without particle rejection (i.e. ? = 1) and those with particle rejection of different rates for the cases of pd = 7 nm and pd = 2 nm are shown in Figure 3.13. These results suggest that the rejection of the copper nanoparticles away from the solid layer into the liquid phase had a 42 negative effect on the movement of the interface, especially for small size particles (Fig. 3.13b). 3.2.2 Temperature distributions in both phases Temperature distributions for both phases will be presented for the case of copper- cyclohexane colloids with initial and surface temperatures of 7.5 oC and 6.0 oC, respectively, with the solid-liquid interface located at s* = 0.002. The effect of decreasing the value of the segregation coefficient (i.e. increasing the particle rejection at the interface) on the temperature distributions in the solid and liquid phases for the case of pd = 7 nm is shown in Figure 3.14. For the case of ? = 0.1, the temperature profiles for the initial volume fraction of particles studied are similar as it is shown in Figure 3.14a. The temperature distributions exhibit a piecewise continuous variation at the interface due to the Stefan thermal condition there (Equation 3.9) and little deviations among the temperature distributions for different initial concentrations of the suspensions is observed. As the value of ? decreased to 0.01, increasing the initial volume fraction of the particles reduces the temperature distributions of both solid and liquid phases on the liquid side of the interface as it is shown in Figure 3.14b. A supercooled layer of liquid next to the interface position (s* = 0.002) extending to x* = 0.003 is clearly observed in Figure 3.14b for an initial volume fraction of 0.01. Reducing the value of ? further (Fig. 3.14c), the temperature distributions are lowered further in the solid and liquid phases for the case of 0? = 0.01. The reduction in temperature is attributed primarily to the reduction of the interface temperature as the value of ? is reduced and the value of 0? is increased. For the suspension with particle size of 2 nm, the temperature distribution especially on the liquid side of the interface are reduced significantly for the higher initial volume fraction values used in this study ( 0? = 0.001 and 0.01), as shown in Figure 3.15. The extent of 43 supercooling in the layer adjacent to the interface is clearly observed in this figure. The effect of particle rejection on the temperature distributions for an initial volume fraction of 0.01 is discussed next. The temperature distributions corresponding to ? = 1 (i.e. no particle rejection) are always higher than those corresponding to lower ? values for both cases of pd = 7 nm and pd = 2 nm shown in Figure 3.16. The temperature in the solid and liquid phase are much lower for ? = 0.001 than ? = 1.0 for pd = 2 nm, as shown in Figure 3.16b. Furthermore, from Figures 16a and 16b, the temperature profiles for ? = 0.01 and ? = 0.001 are nearly identical especially for the pd = 2 nm case, which shows that further reduction in the value of ? may not have a significant effect on the temperature distribution. The dimensionless interface temperature decreases until it reaches an asymptotic value, as the values of ? and pd are reduced (Figure 3.17). This reduction is attributed to the fact that the freezing temperature for the colloid mixture changes with the concentration of the particles. To investigate the existence of the constitutional supercooling during solidification of colloidal suspensions, the variations of the temperature on the liquid side of the interface and the equilibrium temperature are discussed next. For illustration purposes, the instant when the interface is located at s* = 0.002 is considered. For the case of pd = 7 nm and no particle rejection (i.e. ? = 1.0), *LT is always higher than the local equilibrium temperature as shown in Figure 3.18a. This condition will prevail at any location of the interface since the interface temperature for this case is always equal to the initial equilibrium temperature. This behavior is well known to exist for most solidification cases of single-component materials. As the value of ? decreased to 0.1, an extremely narrow region of supercooled liquid next to the accelerated solid-liquid interface (Fig. 3.11a) was created, as observed in Figure 3.18b. As the value of ? 44 was reduced further to 0.01 and 0.001, the region occupied by the supercooled liquid expanded further from the interface as shown in Figures 3.18c and 3.18d, which corresponded to decelerated interface conditions. Figure 3.18a-d shows the region of constitutional supercooled liquid grows with the reduction of the value of ? . As for the conditions with pd = 2 nm, all the cases with particle rejection ( 1?? ) exhibit marked growing zones of constitutional supercooled liquid as shown in Figures 3.19b-d that accompanied decelerated interface movement (Fig. 3.12a-c). In comparing the cases of pd = 7 nm and 2 nm, for the latter case the constitutional supercooled zones have developed much farther away from the interface and the supercooled liquid temperature is lower. The constitutional supercooled area is characterized by the development of dendrite-like structures that have been observed experimentally for colloidal solidification. Furthermore, the zone of the supercooled liquid is one of the factors responsible for the slowdown of the movement of the interface since in all cases where the supercooled zone was developed, a clear slowdown of the interface compared to the pure cyclohexene was observed. 3.2.3 Concentration distributions in both phases Concentration distributions corresponding to copper-cyclohexane with initial and surface temperatures of 7.5 oC and 6.0 oC, respectively, and with the solid-liquid interface located at s*= 0.002 are discussed. The effect of ? on the concentration variations for copper nanoparticles ( pd = 7 nm) is shown in Figure 3.20. Except for the case of no particle rejection ( 1?? corresponding to 1??? ), the concentration profiles follow similar trends. The concentration is constant in the solid phase ( iLiS ?? ???? for **0 sx ?? ), whereas it achieves a maximum value on the liquid side of the interface, and finally the value of the concentration decays until it reaches 45 its initial value. The concentration of the particles at the solid-liquid interface is conserved through Eq. (3.11). It should also be noted that the predicted particle concentration on the liquid side of the interface are many time greater than the initial value. This suggests that the evaluation of the thermophysical properties solely based on the initial values is not a good assumption. For the case of ? = 0.1, the concentration profiles corresponding to different initial volume fractions are similar for all the cases considered, as shown in Figure 3.20a. However, for the cases of ? = 0.01 and ? = 0.001 (Figures 3.20b and 3.20c), the maximum value of the liquid concentration at the interface reduces as the initial volume fraction increases. Lower values of the maximum liquid concentrations correspond to slower solidification rates because the particles had enough time to diffuse away from the interface. It was shown in Figures 3.11b and 3.11c that correspond to the case of pd = 7 nm that the interface slowed down as the initial concentration was increased, which clearly explains the reduction of the peak value of the concentrations of Figures 3.20b and 3.20c. The concentration profiles for the case of pd = 2 nm are shown in Figure 3.21 and the general behavior is similar to that of pd = 7 nm. The value of the maximum concentration on the liquid side of the interface is reduced as the initial concentration is increased. However, for the case of 0? = 0.01, the maximum value of the interface liquid concentration is significantly lower than that of the other initial volume fractions considered. This is due to the extremely slow progression of the solid-liquid interface for this case as it was discussed earlier. This can be attributed to the fact that the value of ? is less than ?t thus promoting control of the solidification process through solutal diffusion, as well as to the constitutional supercooling. 46 The effect of the value of ? on the development of the concentration profiles for the case of 0? = 0.01 is shown in Figures 3.22a-b. As the value of ? decreased, the maximum value of the interface liquid concentration increased, due to the increase of the particles rejected to the liquid side of the interface. For the case of ?pd 2 nm (Fig. 3.22b), particles diffused away from the interface with a less steep slope when compared to the case pd = 7 nm (Fig. 3.22a) and also the maximum value of concentration for the latter case is higher than that of the former case. This can be explained by the higher value of the diffusion coefficient for the case of pd = 2 nm than that of pd = 7 nm, which allows the particles to spread much further from the interface, and thus fewer particles are found ahead of the interface. 3.2.4 Effects of particle properties In order to study the effect of the properties of nanoparticles on the solidification behavior of the NePCM, nanoparticles of three different materials were chosen. The thermophysical properties of these materials are listed in Table 3.1. For the three suspensions of these particles in cyclohexane, the diameter of the nanoparticles was 2 nm. The cold side and initial temperatures were 4.5 oC and 7.5 oC, respectively and 0? = 0.05. Figure 3.23 shows the evolution of the solid-liquid interface for three materials. The movement of the interface is identical for the case of ? = 1 for the three suspensions as it is shown in Figure 3.23a. Extreme, assuming that the rejection rates for the three types of the particles are identical, the value of the segregation coefficient was set to a value of ? = 0.001 to represent a high rejection rate. The development of the interface for the case of ? = 0.001 for the three materials is shown in Figure 3.23b. The interface development is the same for the three suspensions, indicating that for the particles chosen in the current study, the difference in their thermophysical 47 properties did not play a significant role in the solidification process. However, the current model did not take into account the effect of the thermophysical properties on the constitutional supercooling phenomena and the effects of convection were neglected. The variation of the thermophysical properties could play a significant role in these phenomena and influence the solidification process. 3.2.5 Effects of the surface temperature The previous sections showed that the constitutional supercooling becomes more pronounced as the rate of rejection of the nanoparticles is increased, which is responsible for the slowdown of the solid-liquid interface. A simple approach to overcome this slowdown is to increase the value of *T? by decreasing the surface temperature. In this section, the effect of *T? on the constitutional supercooling will be investigated. A copper-cyclohexene suspension was selected for this purpose. The effects of reducing the surface temperature on the movement of the solid-liquid interface and the temperature distribution are presented in Figure 3.24 and Figure 3.25, respectively, for pd = 2 nm and 0? = 0.01. For the range of values of ? , the case of surfT = 2.5 oC exhibits the fastest interface development due to the higher value of *T? . As shown in Figure 3.25, the liquid zone influenced by the constitutional supercooling is extended away from the solid-liquid interface for the case surfT = 2.5 oC if it is compared with the other surface temperature. Two mechanisms compete with each other, supercooling, which helps to speed up solidification, and the constitutional supercooling, which decreases its development. Table 3.3 lists the interface locations at *t = 1000 and 0? = 0.01. The interface is decelerated with a lower rate for the case of surfT = 2.5 oC if it is compared with cases corresponding to surfT = 6.0 oC and 48 surfT = 4.5 oC. 3.2.6 Transition segregation coefficient To illustrate the transition from thermally-controlled solidification to solute-controlled solidification further, consider a specific colloid of cyclohexane and 2 nm diameter copper with Tsurf = 2.5 oC, T0 = 7.5 oC and 0? = 0.01. Based on the thermophysical properties of Tables 1 and 2, the Lewis number is Le= 468. Consider the same colloid except that the solvent is assigned artificially-inflated thermal conductivities in both phases (kS = 2 W/m K and kL = 1 W/m K) that gives rise to a higher Lewis number of Le = 3688. For both colloids, the transitional segregation coefficient remains unchanged at t? = 0.02536. After an elapsed time period, the expediting solidification parameter of the ?enhanced? colloid compared to the actual colloid, as given by '? and defined as: )( )(' nL tyC o n d u c tiv iH ig h e rnL ????? ? (3.28) can be plotted versus the segregation coefficient (Figure 3.26). For values of particle rejection greater than 0.1, the ?enhanced? colloid solidified nearly four times faster than the actual colloid as crystallization is fully controlled by the heat conduction. However, as the value of the segregation coefficient is reduced to 0.01, the expediting factor was lowered markedly, even though the ?enhanced? colloid still solidifies faster than the actual colloid. In this regime, the value of the segregation coefficient is generally higher than t? = 0.02536 and solidification is primarily controlled by diffusion of heat. In effect, due to the inflated thermal conductivity of both phases, the ?enhanced? colloid crystallized faster than actual colloid. However, as the value 49 of the segregation coefficient was reduced below the value of the transitional segregation coefficient, '? attained an asymptotic value of nearly 1.25, which indicates that the two colloids freeze nearly at the same rate. In this regime, solidification is controlled by the diffusion of the solute and the artificial improvement of thermal conductivity has no significant influence on the crystallization process. One can conclude from Figure 3.26 that the main influential parameter for harnessing the desired expedited freezing of NePCM is the transitional segregation coefficient rather the effective thermal conductivity. This quantity is the governing parameter that dictates if the benefits of the high effective thermal conductivity reported in [Yavari et al. (2011), and Zheng et al. (2011)] can be utilized in colloid-based thermal energy storage applications that involve latent heat. 3.3 Comparison with the Experimental Results The predicted freezing times for the interface to travel a certain distance obtained for the one-dimensional model were compared with the experimental data of Fan and Khodadadi (2011) for case of cyclohexene-copper oxide suspensions, as shown in Figure 3.27. The cold surface temperature is maintained at -15 oC. The diameter of the particles is chosen to be 7 nm, which is very close to the value reported by Fan and Khodadadi (2011) of 8 nm. The results are obtained for different initial volume fractions, and three different values of the segregation coefficient. Knowing the exact value of the segregation coefficient for each case is impossible since there is no theoretical formulation for the segregation as a function of the operational parameters. Experimental observations indicate that as the volume fraction of the particles increases, the amount of the particles rejected out of the crystalized phase is increased, and thus we anticipate that the predicted results of the lower segregation coefficients will match those of high-volume fractions. As shown from Figure 3.27, the freezing times obtained from the analytical solution 50 and the experiments are close for the pure cyclohexene case. The results also show that as the volume fraction increases, the experimental and analytical results are closer as the value of the segregation coefficient is decreased. 3.4. Summary By using the closed-form analytical solution of an extended Rubinstein problem, an analysis of the effect of the rejection of the particles on the solidification process of a NePCM was conducted, and the following conclusions are drawn: 1- For the case of no particle rejection, it should not be assumed that through addition of more nanoparticles to increase the thermal conductivity, expedited freezing compared to the pure PCM will occur. This is due to the behavior of the dimensionless temperature difference of the initial field and the cold surface which exhibits inverse dependence on both particle size and concentration. 2- When particles were rejected from the crystalline phase, two mechanisms are controlling the solidification process. One was thermal diffusion, and the other was the diffusion of the solute. The transition segregation was the controlling parameter that dictates the importance of each mechanism during the process. 3- For the cases of ? < 1.0 (i.e. particle rejection) when ? > ?t expedited freezing with respect to the pure PCM occurs. Otherwise, the interface will decelerate because the solutal mechanism would overtake the thermal mechanism. 4- For the cases when ? < t? the interface was decelerated with respect to the pure PCM as the volume of the particles was increased, and this deceleration is more pronounced as the particle size was reduced. 51 5- As the value of the segregation coefficient was reduced the maximum concentration at the solid-liquid interface is reduced as the concentration of the particles is increased, and as the particle size was reduced. 6- As the value of the segregation coefficient was reduced a zone of constitutional supercooling wasformed in front of the solid-liquid interface, and this zone grew as the particle size was reduced. 7- The transition segregation coefficient was decreased as the surface temperature was reduced, and it also increased as the volume fraction of the particles is increased, and the particle size was educed. 52 Table 3.1: Thermophysical properties of cyclohexane Liquid Solid Specific Heat (J/kg K) 1762.8 1800 Density (kg/m3) 779 856 Thermal Conductivity (W/m K) 0.127 0.1359 Dynamic Viscosity (Pa s) 0.974x10-3 - Latent Heat (J/kg) 32,557 Melting Temperature (K) 279.5 53 1Table 3.2: Thermophysical properties of the nanoparticles Al2O3 CuO Cu Density (kg/m3) 3880 6510 8954 Specific Heat (J/kg K) 729 540 383 Thermal Conductivity (W/m K) 42.34 18 400 54 Table 3.3: Solid-liquid interface locations for *t = 1000, 0? = 0.01 and different values of ? * 1=1000,=* ?ts * 0.001=1000,=* ?ts 100)( * 1=1000,=* * 0 . 0 0 1=1000,=** 1=1000,=* ?? ? ?? t tt s ssA b s Tsurf = 6.0oC 0.0016 0.0007 56.25% Tsurf = 4.5oC 0.0039 0.0019 51.25% Tsurf = 2.5oC 0.0057 0.0028 50.87% 55 Figure 3.1 The Le and St numbers with 0? and pd for cyclohexane-based colloids of (a, b) CuO, (c, d) Al2O3 and (e, f) Cu. (a) (b) (c) (d) (e) (f) 56 (a) (b) (c) Figure 3.2 Variation of *T? with 0? and pd for cyclohexane-based colloids of (a) CuO, (b) Al2O3 and (c) Cu 57 Figure 3.3 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 7 nm) for various values of 0? (no particle rejection at the interface) 58 Figure 3.4 Variation of the solid-liquid interface location with time for the case of copper nanaoparticles ( pd = 5 nm) for various values of 0? (no particle rejection at the interface) 59 Figure 3.5 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 2 nm) for various values of 0? (no particle rejection at the interface) 60 Figure 3.6 Comparison among the variations of the solid-liquid interface with time for different copper nanoparticle sizes and 0? = 0.01 (no particle rejection at the interface) 61 Figure 3.7 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 2 nm) and various values of 0? for 0T = 7.5oC and surfT = 2.5oC (no particle rejection at the interface) 62 Figure 3.8 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 2 nm) and various values of 0? for 0T = 7.5oC and surfT = 0.5oC (no particle rejection at the interface) 63 Figure 3.9 Variation of the nondimensional supercooling with the 0? for different values of surfT and pd = 2 nm 64 Figure 3.10 Variation of ? with surfT for different values of 0? and pd 65 Figure 3.11 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 7 nm) and various values of 0? (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001 (a) (b) (c) 66 Figure 3.12 Variation of the solid-liquid interface location with time for the case of copper nanoparticles ( pd = 2 nm) and various values of 0? (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001 (a) (b) (c) 67 (a) (b) Figure 3.13 Variation of the solid-liquid interface location with time for the case of copper nanoparticles with 0? = 0.01 and different rejection rates for (a) pd = 7 nm and (b) pd = 2 nm 68 Figure 3.14 Variation of temperature with *x for the case of copper nanoparticles ( pd = 7 nm) with *s = 0.002 and various values of 0? for (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001 (a) (b) (c) 69 Figure 3.15 Variation of temperature with *x for the case of copper nanoparticles ( pd = 2 nm) with *s = 0.002 and various values of 0? for (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001 (a) (b) (c) 70 (a) (b) Figure 3.16 Variation of temperature with *x when *s = 0.002 for different ? values with 0? = 0.01 for (a) pd = 7 nm and (b) pd = 2 nm 71 Figure 3.17 Variation of the interface temperature iT* with ? for different values of pd 72 (a) (b) (c) (d) Figure 3.18 Variation of *LT and *equT with *x for the case of pd = 7 nm when *s = 0.002 and 0? = 0.01 for (a) ? = 1.0, (b) ? = 0.1, (c) ? = 0.01 and (d) ? = 0.001. 73 (a) (b) (c) (d) Figure 3.19 Variation of *LT and *equT with *x for the case of pd = 2 nm when *s = 0.002, 0? = 0.01 for (a) ? = 1.0, (b) ? = 0.1, (c) ? = 0.01 and (d) ? = 0.001. 74 Figure 3.20 Variation of concentration with *x for copper nanoparticles ( pd = 7 nm) with *s = 0.002 and various values of 0? for (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001. (a) (b) (c) 75 Figure 3.21 Variation of concentration with *x for copper nanoparticles ( pd = 2 nm) with *s = 0.002 and various values of 0? for (a) ? = 0.1, (b) ? = 0.01 and (c) ? = 0.001. (a) (b) (c) 76 (a) (b) Figure 3.22 Variation of concentration with *x with *s = 0.002 for different ? values for the case of 0? = 0.01: (a) pd = 7 nm and (b) pd = 2 nm 77 (a) (b) Figure 3.23 Variation of the solid-liquid interface location with time for three materials for the case of 0? = 0.05 and pd = 2 nm: (a) ? = 1.0 and (b) ? = 0.001. 78 (a) (b) (c) (d) Figure 3.24 Variation of the solid-liquid interface location with time for the case of 0? = 0.01, pd = 2 nm and for various *surfT values: (a) ? = 1.0, (b) ? = 0.1, (c) ? = 0.01 and (d) ? = 0.001. 79 (a) (b) (c) (d) Figure 3.25 Variation of temperature with *x for *s = 0.002 for the case of 0? = 0.01, pd = 2 nm and for various *surfT values: (a) ? = 1.0, (b) ? = 0.1, (c) ? = 0.01 and (d) ? = 0.001. 80 Figure 3.26 Variation of the expedited freezing parameter epsilon(?) with the segregation coefficient (?). 81 Figure 3.27 Comparison between the predicted freezing times based on the 1-D model and experimental data of Fan and Khodadadi (2011). 82 Chapter 4 Two-Dimensional Model for Colloidal Solidification In this chapter, the one-dimensional Rubinstein problem for solidification of colloid suspensions has been extended to two dimensions, in effect leading to fluid flow. For this purpose, a two-dimensional one-fluid mixture model combined with the enthalpy-porosity model has been used. The model takes into account both the Brownian diffusion and thermophoretic effects. Using a water-copper NePCM as the working colloid, and with the diameter of copper nanoparticles having two distinct values of 5 nm and 2 nm, the NePCM will be solidified from the bottom side of a cavity unidirectionally. Results of the development of the liquid fraction, concentration and the liquidus and solidus temperatures will be presented and discussed. 4.1 Mathematical Modeling Since colloidal suspensions can be considered as binary mixtures, their solidification processes are accompanied by formation of a mushy zone. A mushy zone can be defined as a region where both liquid and solid phases coexist and the phase change process occurs within its boundary. The model that will be used for this study is based on the one-fluid mixture model that has been used for modeling the binary alloy solidification, by Bennon and Incropera (1987), Beckermann and Viskanta (1988), and Voller and Prakash (1987). 4.1.1 Governing equations The volume-averaged equations are obtained through integration of the transport equations over a very small control volume. The volume-averaged transport equations are valid 83 in the liquid and solid phases as well as in the mushy region. The following assumptions have been used for deriving the volume-averaged transport equations: 1. Two-dimensional Cartesian coordinate system, 2. The base material in the liquid and solid phases is homogenous/isotropic and their properties are identical, 3. The solid and liquid in the mushy zone are in local equilibrium, 4. The solid phase is stationary, 5. The Boussinesq approximation is used to account for buoyancy-driven thermal-solutal convection effects due to the variation of the density with temperature and concentration; however, due to the thermally-stable nature of the stratified system studied, only solutal convection will be considered. 6. The thermophysical and transport properties of the dilute ( 1??? ) colloid suspension are not constant and change with the concentration of the nanoparticles, 7. Velocities due to density change upon phase change are neglected. 8. The suspended particles are assumed to be mono size and agglomeration of the particles subsequent to initiation of freezing is ignored. The pertinent volume-averaged dimensional transport equations are: Continuity: ,0U??? ? (4.1) 84 Momentum: ,)()()1()( 3 22 yr e fm eTTgUCUpUUtU ?????? ?????????????? ??? ????? (4.2) Thermal Energy: ,tL)Tk(TUctTc pp ???????????? ???? ? (4.3) Species: ),T TDD(Ut TwBww ??????????? ??? ? (4.4) with the local values of the particle weight fraction ( w? ) and colloid liquid fraction (? ) related through the following relation: .)1( wswlw ????? ??? (4.5) The momentum equation contains the Darcy Law damping terms and the associated porosity constant (Cm) was set to 105 kg/m3s. Moreover, the contribution of nanoparticle diffusion to the heat flux vector in the thermal energy equation is neglected (Buongiorno, 2006). The governing equations are non-dimensionalized using the following non-dimensional groups: T T B B CHp m CHp CH C D Le D Le L TTc St e C Da HTTg Ra k c TT TT H Y y H X x H tHU u ?? ? ?? ? ? ? ?? ? ? ? ?? ? ? ? ? ? ?? ? ? ?????? ,, )( , )1( , )( )( ,Pr ) ( ( ,,,, 33 2 ? ? (4.6) 85 The non-dimensional volume- averaged equations are listed below: Continuity: ,0??? u? (4.7) Momentum: yeRaDa uupuuu ? ????? PrPrPr)( 2 ????????????? ? (4.8) Thermal Energy: ??? ??????????????? S teu 1)(? (4.9) Species: ]11[)( ????????????? TwBw w LeLeu ???? ? (4.10) 4.1.2 Initial/Boundary conditions The geometry of the physical model considered in the current investigation is shown in Figure 4.1 with the origin of the adopted Cartesian coordinates at the lower left corner. The motionless colloid is contained in a square cavity (side H) with an initial temperature higher than that of the liquidus temperature corresponding to an initial concentration. Thus: .0tf o r1.0,K274T,0U w i ninin ???? ?? (4.11) For 0t? , the dimensional boundary conditions for the current problem are: 86 .Hy0,Hxf or,0 x ,0 x T ,0U ,Hy0,0xf or,0 x ,0 x T ,0U ,Hy,Hx0f or,0 y ,K273T,0U ,0y,Hx0f or,0 y ,K268T,0U w w w w ???? ? ? ? ? ? ? ???? ? ? ? ? ? ? ???? ? ? ?? ???? ? ? ?? ? ? ? ? ? ? ? ? (4.12) Thus, solidification of the colloid began by lowering the temperature of the bottom wall below that of the liquidus temperature and the liquid/solid interface advanced into the liquid colloid. 4.1.3 Mixture relations and effective thermophysical/transport properties It is assumed that the particles suspended in the colloid are hard spheres and no clustering takes places. The properties of the base fluid are assumed to be identical for the solid and liquid phases. The thermophysical and transport properties of the colloid vary with the volume fraction (? ) of the particles. The colloid density, heat capacities and part of the Boussinesq terms are calculated from the mixture relation: ,)1( pf ????? ??? (4.13a) ,)c()c)(1(c ppfpp ????? ??? (4.13b) .)())(1( pf ???????? ??? (4.13c) The mass fraction can be converted to the corresponding volume fraction via the following relation: 87 .)1( wpwf wf ???? ??? ??? (4.14) The viscosity of the colloid is obtained from: .)1( 5.2f??? ?? (4.15) The thermal conductivity of the colloid suspension is evaluated using: ,21 kkk ?? (4.16a) That combines the Maxwell relation: ,)(2 )(221 pffp pffpf kkkk kkkkkk ??? ???? ?? (4.16b) and the enhancement due to thermal dispersion given by: .dU)c(Ck ppk2 ?? ?? (4.16c) The quantity 51.0?kC stands for an empirical constant obtained from Wakao and Kaguei (1982). To calculate the mass diffusivity (Brownian diffusivity) of the nanoparticles in the colloid suspension, the compressibility factor )(z? that accounts for the effect of particle-particle 88 interaction on the osmotic pressure was calculated from the following empirical fitting relation obtained from Peppin et al.(2006): . 64.01 )64.0 1018()64.0 410()64.0 14(1 )( 32 ? ??? ? ? ?????? ?z (4.17a) The mass diffusivity of the particles that is a function of the volume fraction of the nanoparticles is then calculated from: ,d )z(d)1(DD 60B ????? (4.17b) Where 0D is calculated from the Einstein-Stokes relation: .d3TkD p B0 ??? (4.17c) The thermophoretic diffusivity is computed from the following relation: ,D kT ????? (4.18) Where pf fk kk2 k26.0 ??? is a non-dimensional parameter that is a function of the thermal conductivities of base liquid and the nanoparticle as it is suggested by Buongiorno [2006]. Sensible enthalpy and latent heat of fusion make up the enthalpy of the colloid and the latent heat is calculated from the following relation: .)L)(1()L( f??? ?? (4.19) 89 The liquidus and solidus temperatures of the colloidal system are calculated from Equations (4.20a) and (4.20b) that are linear models of the typical temperature vs. concentration relation of phase diagram of a binary mixture: ,wlmLiq mTT ??? (4.20a) , 0 w lmSol kmTT ??? (4.20b) Where Tm, ml and k0 represent the melting temperature of the pure solvent (in this case water taken as 273 K), the liquidus slope and the segregation coefficient, respectively. The segregation coefficient represents the ratio between mass fraction of the particles on the solid and liquid sides of the solid-liquid interface. The liquidus slope is given by the following relation that was developed specifically for colloidal suspensions as given in Peppin et al. (2006): .Lv Tkm ffp 2mB l ?? (4.21) A value of 0.1 was assigned to the segregation coefficient (k0) that was kept constant throughout the simulation. This can be justified due to the small difference between the hot and cold wall temperatures that promotes a slow solidification rate. This, in turn, ensures that equilibrium conditions prevail and a linearized model for the phase diagram is valid. The colloid liquid fraction (? ) is related to the local liquidus and solidus temperatures via the following relations: 90 .1 , ,0 L i q L i qSol SolL i q Sol Sol TTfo r TTTfo rTT TT TTfo r ?? ????? ?? ? ? ? (4.22) Similar to Beckermann and Viskanta (1988) and adopting ?,H and (TH ? TC) for forming dimensionless groupings, the Brownian Lewis (LeB), Thermophoretic Lewis (LeT), Darcy (Da), Prandtl (Pr), Rayleigh (Ra) and Stefan (Ste) numbers govern this problem (see relations 4.8-4.10). The mixture model utilized here has been extensively used for simulating binary alloy solidification problems [Bennon and Incropera (1987), Beckemann and Viskanta (1987)]. In the present study, the applicability and relevance of this model is extended to solidification of colloidal suspensions (Peppin et al. (2006)) of nano-scale particles. Specifically, the effects of the diameter of the spherical particles (Equations 4.16c, 4.17c, 4.21) on modifying the concentration, liquidus and solidus temperature fields will be assessed. 4.2 Solution Procedure The finite volume method was used to discretize the governing equations. The grid was staggered and distributed uniformly in the x and y directions. The SIMPLE algorithm was used to handle the coupling of the pressure and the velocity field. The second-order upwind scheme was used to handle the coupling of convection and diffusion terms in the momentum, energy and species equations. The enthalpy porosity method was used to simulate the solidification process. The implemented convergence criterion required that at the end of each time step the residuals of the momentum equation be lower than 10-5, whereas for energy and species equations this value 91 was 10-7. The variable time step was set to 0.01 s for 1000 time steps, then increased to 0.05 s for another 1000 time steps and finally raised to 0.1 s for the remainder of the simulation. The computational fluid dynamics software FLUENT 6.3 was used to implement the details of the numerical model discussed above. Special user-defined functions were developed to handle the variation of the thermophysical and transport properties with the volume fraction of the nanoparticles. A grid independence study was conducted using the grid arrangements of 50 x 50, 100 x 100 and 150 x150. The variations of the temperature, liquid fraction ( ) and mass concentration of the nanoparticles ( ) at x = 0.5 H at two times (t = 10 and 100 s) for the case of dp = 5 nm particles are shown in Figure 4.2. At each time, the variations of the monitored quantities corresponding to different grid systems are banded and the prevailing trends are recovered. However, the 100 x100 grid system was selected due to its accuracy and reduced simulation time. The results of the adopted model were also compared to those of Hannoun et al. (2005) for the case of melting of pure tin inside a rectangular cavity heated on the vertical walls. Comparing Figures 4.3(a) to 4.3(b) and 4.3(c) to 4.3(d), it is observed that the present model adequately captures the presence of a single recirculating cell at time instant t = 100 s and its breakdown to a multi-cell structure at t = 200 s. More important, the tabulated data on the positions of the liquid-solid interface at the same time instants Hannoun et al. (2005) matched very closely to the current predictions (Figure 4.3e). 4.3 Results and Discussion The numerical investigations were carried out for the case of unidirectional solidification of water with copper nanoparticles within a H = 10 mm square. The initial temperature of the colloid was 274 K with a 10% nanoparticle mass fraction (1.22 vol%) and the reported ? w? 92 nanoparticle diameters were 2 and 5 nm. Other particle diameters in the 2-10 nm range were also explored since NePCM freezing tests for such colloids have been reported by Fan and Khodadadi (2012). Adhering to the utilized binary alloy model that considers a colloid as a solution is another reason for focusing on nanoparticles with single-digit diameters. The thermophysical properties of the solvent and the nanoparticles are listed in Table 4.1 and the corresponding values of the pertinent dimensionless groupings at t = 0 are summarized in Table 4.2. Furthermore, the values of the liquidus slope and segregation coefficient are listed in Table 4.3. The instantaneous liquid fraction field (? ) and the transient development of the solid- liquid interface for the case of dp = 5 nm are shown in Figure 4.4. Initially, the liquid colloid (red zone) completely occupies the domain of computation. At early stages of the solidification process observed in Figures 4.4(a) and (b), the advancing interface has a planar shape and the thickness of the mushy zone is small. At later time instants (Figures 4.4(c) and (d)), the advancing liquid-solid interface still maintains its nominal planar shape and the thickness of the mushy zone has increased slightly. This indicates that during the solidification process, the difference between the solidus and liquidus temperatures was small which, can be justified by the low value of the liquidus slope for the current case (Table 4.3). Moreover, the nanoparticles did not redistribute greatly on the liquid side of the interface at early stage of solidification, as shown below. The development of the solid-liquid interface and the instantaneous liquid fraction field for the case of smaller nanoparticles (dp = 2 nm) are shown in Figure 4.5. Early on at t = 10 s, the interface is planar and is similar to the previous case (Fig. 4.4a). However, as time proceeds to t = 100 s, appearance of repeating cell-like pockets at the interface is observed. Furthermore, the thickness of the mushy zone increases. This can be attributed to the rejection of 93 the nanoparticles away from the thickening solid layer into the shrinking liquid layer. In turn, the increased concentration of the particles ahead of the mushy zone leads to regions of different liquidus and solidus temperatures within the mushy zone that promotes faster crystallization of parts of the suspension. This is because the liquidus slope had a greater value than the previous case (0.735 compared to 0.046 K) as dictated by the inverse cubic dependence of Equation (4.15) on the particle diameter. According to Equations (4.14a-b), the liquidus temperature is more responsive to the concentration change of the nanoparticles. As time proceeds as in Figures 4.5(c) and (d), the repeating cellular pockets grow in size to take the shape of a dendritic structure of the order of 0.1 H. Similar shapes have been observed experimentally for the case of alumina-water colloidal suspension (Deville et al., 2009). Simulations of the same problem with horizontal cavity lengths of 0.5 H and 1.5 H were also performed so as to rule out any possible dependence of the observed structures on artificial computational inaccuracies. The mechanism that gave rise to the above-mentioned dendrite-like structures is generally referred to as constitutional supercooling. Constitutional supercooling ensues upon lowering of the liquidus temperature ahead of the interface due to the increased particle concentration in that area due to the particle rejection and its redistribution in the liquid phase. In the following section, the evolving mass concentration field, liquidus and solidus temperature distributions are discussed in an attempt to explore the mechanisms behind the development of dendritic or non-dendritic interfaces during the solidification of colloidal suspensions. Colorized contours of the mass concentration field of the nanoparticles ( w? ) during the solidification process of the water-copper colloid for the case of dp = 5 nm are given in Figure 4.6. It is clearly observed from Figures 4.6(a) and (b) that the particles were rejected from the growing solid phase and have accumulated on the liquid side of the planar-like interface creating 94 a thin layer of highly concentrated nanoparticles. As time proceeds as in Figures 4.6(c) and (d), the concentrated layer of the nanoparticles created by the rejected particles becomes thicker, leading to discrete pockets along the interface where the concentration field achieves its local maximum values as much as 4.5 times that of the initial concentration (red zones). However, the rejected and redistributed nanoparticles did not diffuse markedly away from the interface. Concurrent to formation of high concentration pockets of particles, regions of low nanoparticle concentration with varying volumes are left behind within the solid layer below the interface shown as dark blue regions in the figure. The development of the mass concentration field of the nanoparticles for the case of dp = 2 nm is shown in Figure 4.7. At time t = 10 s, the rejected particles form a layer above the thickening solid phase. However, as the time proceeds to t = 100 s, the layer of the segregated particles becomes thicker if it is compared with the t = 10 s snapshot (Fig. 4.7a), as well as at the same time for the case of dp = 5 nm (Fig. 4.6b). This indicates that the Brownian diffusion of nanoparticles is more important than in the previous case. At t = 500 s, the nearly uniform layer of the nanoparticles has transformed into discrete cells of high concentrations of nanoparticles, separated by regions of low concentrations. Similar observations were reported by Deville et al. (2009) for the case of the experimental solidification of colloidal suspensions. At t = 1000 s, cellular patterns of regions of high concentration are clearly observed. What distinguishes Fig 4.7(d) is that nanoparticles diffuse away from the interface similar to a rising plume that creates a concentration gradient in the mushy and liquid zones. This concentration gradient is responsible for the genesis of the dendrite structures that were discussed above in Figure 4.5. To provide greater insight into the complex transport of the particles in the vicinity of the liquid-solid interface, an isometric view of the mass concentration field corresponding to the case of dp = 2 nm at t = 1000 s is shown in Figure 4.8. At the given 95 time instant, the uniformly distributed nanoparticle concentration field is found at the initial value of w? = 0.1 for 5.0/ ?Hy . Zones of depleted particle concentration ( 3.0/ ?Hy ) that are left behind the advancing interface along with somehow regularly spaced, nearly identical regions ahead of the interface that have been enriched by the rejected particles are clearly observed. Based on the assumptions of the model utilized in this study, there are two mechanisms that are responsible for the redistribution of the nanoparticles. These are the Brownian diffusion and convection generated from thermal-solutal effects. A close-up view of the liquid fraction field in the vicinity of the interface that is superimposed with the predicted fluid velocity vectors for case of dp = 2 nm at t = 1000 s is shown in Figure 4.9. In effect, a small portion of the interface ( 8.0/3.0 ?? Hx ) from Figure 4.5 has been reproduced. Only six of the dendritic structures are observed in the close-up view and the convection cells do not extend beyond the interface zone. A multitude of recirculating fluid cells (solutal convection) due to the density gradient in the vicinity of the interface are observed. The maximum observed velocities for the case of dp = 5 nm and dp = 2 nm were 3.82x10-5 and 6.72x10-5 m/s, respectively. Concurrent to solutal convection, the particle diffusion mechanism that is also responsible for the movement of the nanoparticles manifests itself through the Brownian diffusion coefficient that depends inversely on the particles size (Equation 4.11c). For the two cases studied here, the Brownian diffusivity was approximately doubled when the size was reduced from 5 nm to 2 nm and the corresponding Brownian Lewis numbers were 1,660 and 664, respectively (Table4.2). Variations of the liquidus temperature fields during solidification are illustrated in Figure 4.10 for the case of dp = 5 nm. A layer of low value of the liquidus temperature is created above the crystallized portion of the suspension as it is shown in Figures 4.10(a) and (b) due to the high concentration of the nanoparticles in that region. At later times (shown in Figures 4.10(c) and 96 (d)), the maximum and minimum liquidus temperatures at t = 1000 s are 272.979 and 272.998 K, respectively. This indicates the liquidus temperature is fairly uniform throughout the suspension and almost identical to the freezing temperature of pure water (Tm = 273 K). This was attributed to the small value of the liquidus slope (i.e. -0.047 K) that reflects a weak relation between the concentration of the nanoparticles and the liquidus temperature. Similar trends can be seen for the case of the development of the solidus temperature for the case of dp = 5 nm that is shown in Figure 4.11. A layer of lower solidus temperature is observed above the solid-liquid interface. The maximum and minimum values of solidus temperature at t = 1000 s are 272.792 and 272.976 K, respectively. These values are close to each other and are nearly equal to that of the melting temperature of the pure solvent (water for the current case). Since the values of the liquidus and solidus temperatures did not vary significantly throughout the suspension, this indicates that the suspension will be solidified uniformly at the same temperature, which is one of the reasons that dendritic structures did not appear for the dp = 5 nm case. The instantaneous variations of the liquidus temperature for the case of dp = 2 nm are shown in Figure 4.12. For t = 10 and t = 100 s, a layer of low value of the liquidus temperature has emerged above the solid phase due to the high numbers of rejected nanoparticles in this region. At the same time, a layer of high liquidus temperature was created below the interface due to the lower nanoparticles content there. At later time instances, i.e. t = 500 s and t = 1000 s, one can observe distinct zones of high and low liquidus temperatures. For example, the maximum and the minimum values of the liquidus temperature for the case of t = 1000 s are 97 272.81 and 272.979 K, respectively. The difference between these two values is about 0.169 oC, which is about 9 times greater than the value (0.019 oC) discussed above in Figure 4.10(d) for dp = 5 nm. At the early stages of the solidification process, the evolution of the solidus temperature was similar to that of the dp = 5 nm. As shown in figure 4.13. However, at later time instants, distinct areas of high and low values of the solidus temperature evolved. The difference between the minimum and maximum solidus temperatures at those time instances for the case of dp = 2 nm is greater than in the previous case (dp = 5 nm). For example, for the case of t = 1000 s, the minimum and maximum solidus temperatures are 271.1 and 272.7 K, respectively. The difference being about 1.6 oC is greater that the dp = 5 nm case that was equal to 0.184 oC. This noticeable difference in the solidus temperatures is a driving force for parts of the suspension to crystallize faster than the other portions, thus leading to creation of dendrites. In Figure 4.14, the developing profiles of the nanoparticle mass concentration at x = H/2 are shown for the case of dp = 5 nm at different time instants. As the interface sweeps through the colloid, the particle concentration rises above the initial value of 0.1 ahead of the interface due to the rejection of the nanoparticles from the solid phase and the concentration of the particles in the solidified layer drops below the initial value. The particle concentration then reaches a maximum value and decreases until it attains the initial concentration value of 0.1. At the monitored position (x = H/2), values of the highest observed concentration increase at later time. This is clearly shown by the observed peaks on the graph that indicates that the nanoparticles have accumulated ahead of the interface. 98 The instantaneous nanoparticle mass concentration profiles for the case of dp = 2 nm on the vertical mid-plane (x = H/2) and both sides of it are shown in Figure 4.15. Different vertical planes are chosen along a growing dendrite as well as on its left and right sides. The positions of these vertical lines are superimposed on the mass concentration field at t = 1000 s (the same as Figure 4.7d) at the top of Figure 4.15. The transient evolution of the nanoparticle mass concentration along the left side of a dendrite is shown in Figure 4.15a. At early times t = 10 and t = 100 s, the concentration profiles are similar to that of the dp = 5 nm case, however, the values of the maximum concentrations are smaller than those of the dp = 5 nm case (Fig. 4.14). At later times, t = 500 and t = 1000 s, the concentration profiles exhibit both global and local maxima. This may be attributed to the evolving solutal convection effects. The particle concentration along a dendrite on the vertical mid-plane (x = H/2) is shown in Figure 4.15b. At t = 10 and 100 s instants, the trends of the concentration profiles are similar to that for the dp = 5 nm, exhibiting rise to a maximum value and then dropping to the initial value of 0.1. However, at the time instants t = 500 and t = 1000 s, the concentration for a significant part of the cavity height is below the initial mass concentration (0.1) that indicates a zone of significant particle depletion/rejection. The location of the profile of Figure 4.15b was specifically chosen along a growing dendrite to illustrate clearly that the dendrite is the site in the growing solid where the particle concentration is low compared to the other locations in the suspension. The other intriguing feature is that the nanoparticles do not accumulate in front of the dendrite tip and this can be explained from the low values of the maximum concentrations at times t = 500 and 1000 s. The concentration profiles on the right side of the dendrite (Figure 4.15c) are similar to those of Figure 4.15a. However, at t = 1000 s, the global and the local maxima are observed more clearly, which may be an indication that solutal convection is more active in this region than on 99 the left side of the mid-plane. Furthermore, comparison of Figures 4 14 and 4.15b suggests that the growing solid along the vertical line x=H/2 for the case of dp = 2 nm consists of fewer nanoparticles. This leads one to conclude that the solid phase resulting from a dendritic growth will consist of zones of less concentration of nanoparticles compared to one resulting from a planar interface at the same operating conditions and even with equal segregation coefficients. In relation to conservation of particles over the solution domain, it should be noted that the concentration field was integrated at every time instant and the deviations from the initial concentration field was determined to be negligible. The instantaneous liquid fraction fields at 1000 s for pure water (i.e. w? = 0) and those of w? = 10% colloids with 5 and 2 nm diameter nanoparticles are given in Figure 4.16. The thickness of the frozen layers for pure water and w? = 10% (dp = 5 nm) colloid were nearly identical. It is reasoned that the expedited freezing accommodated by an increase in the thermal conductivity achieved by adding of the nanoparticles (as predicted by a thermal model limited to static particles (Fan and Khodadadi, 2012) was balanced by a reduction in the liquidus and solidus temperatures. However, for the case of dp = 2 nm nanoparticles after 1000 s, the frozen layer is thinner when compared to that of pure water. This can be attributed to nanoparticle redistribution that is responsible for lowering the solidus temperature as much as 2 oC below that of pure water. Even though the thermal conductivity was increased in the particle-rich pockets, the solidification rate was not expedited. One can conclude that for binary mixtures or colloidal solidification, increasing the thermal conductivity does not guarantee faster solidification rates. 100 4.4 Summary A numerical investigation of the solidification process of a water/copper nanoparticles colloid model that considered the suspension as a binary mixture was conducted and the following conclusions were drawn: 1. While the interface was planar for the 5 nm nanoparticles , as the particle size decreased, the interface changed from a planar stable morphology to an unstable dendritic structure. The transition was linked to the widening of the difference between the minimum and maximum values of the liquidus or solidus temperature fields. 2. Redistribution of the rejected nanoparticles from the liquid-solid interface by the Brownian diffusion and solutal convection were observed to be more marked for the smaller particle size of 2 nm. 3. For the case of 5 nm nanoparticles, the rejected particles accumulated in front of the interface, whereas for the case of 2 nm nanoparticles, dendritic structures were created and the particles accumulated between these cells. This was attributed to more noticeable effect of the constitutional supercooling for smaller size nanoparticles whereby the solidus temperature was lowered. 4. The solid structure that evolved from a dendritic growth contained zones of lower nanoparticle concentration than the other parts of the suspension due to size-dependent rejection of nanoparticles. 5. By decreasing the size of the nanoparticles, the layer of frozen colloid that formed in a given time was thinner suggesting that expedited freezing due to higher thermal conductivity was negated by constitutional supercooling. 101 4.1 Thermophysical and transport properties for the solvent and the nanoparticles Water Copper nanoparticles Density 997.1 kg/m3 8954 kg/m3 Viscosity 8.9x10-4 Pa s - Specific Heat 4179 J/kg K 383 J/kg K Thermal Conductivity 0.6 W/m K 400 W/m K Thermal Expansion Coefficient [What Ref.? 4] 2.1x10-4 K-1 1.67x10-5 K-1 Heat of Fusion 3.35x105 J/kg - 102 Table 4.2 Values of the dimensionless groupings at the onset of solidification (t = 0) Brownian Lewis Number BD ? Thermophoretic Lewis Number TD ? Prandtl Number kcp? Rayleigh Number ???? 3)()( HTTg CH ? Stefan Number L )TT(c CHp ? w? = 0 - - 6.2 8x104 0.623 dp =2 nm 664 3.75x104 5.6 7.44x104 0.63 dp =5 nm 1,660 3.75x104 5.6 7.44x104 0.63 103 Table 4.3 Liquidus slopes and segregation coefficient used in the current study dp = 5 nm dp = 2 nm Liquidus slope 0.047 K 0.735 K Segregation Coefficient 0.1 0.1 104 Figure 4.1 Geometry of the physical model 105 ` Figure 4.2 Grid independence results showing (a) temperature, (b) liquid fraction and (c) mass fraction of the nanoparticles at x = 0.5 H for different grids (dp = 5 nm) at two time instants. 106 Figure 4.3 Comparison between the results of the current model and Hannoun et al. (2005): (a) stream function at t = 100 s Hannoun et al. (2005), (b) stream function at t = 100 s (current model), (c) stream function at t = 200 s Hannoun et al. (2005), (d) stream function for t = 200 s (current model) and (e) instantaneous positions of the liquid-solid interfaces (a) (b) (c) (d) (e) 107 (a) (b) (c) (d) Figure 4.4 Development of the liquid fraction field (? ) at different time instants: (a) 10 s, (b) 100 s, (c) 500 s and (d) 1,000 s for an initial mass concentration field of 10% and dp = 5 nm. 108 (a) (b) (c) (d) Figure 4.5 Development of the liquid fraction field (? ) at different time instants: (a) 10 s, (b) 100 s, (c) 500 s and (d) 1,000 s for an initial mass concentration field of 10% and dp = 2 nm. 109 (a) (b) (c) (d) Figure 4.6 Development of the nanoparticle concentration field ( w? ) at different time instants: (a) 10 s, (b) 100 s, (c) 500 s and (d) 1,000 s for an initial mass concentration field of 10% and dp = 5 nm. 110 (a) (b) (c) (d) Figure 4.7 Development of the nanoparticle concentration field ( w? ) at different time instants: (a) 10 s, (b) 100 s, (c) 500 s and (d) 1,000 s for an initial mass concentration field of 10% and dp = 2 nm 111 Figure 4.8 Isometric view of the nanoparticle concentration field ( w? ) at t =1,000 s for an initial mass concentration field of 10% and dp = 2 nm. ?w 112 Figure 4.9 Velocity vectors in the vicinity of the liquid-solid interface ( 8.0/3.0 ?? Hx ) at 1,000 s for an initial mass concentration field of 10% and dp = 2 nm. 113 Figure 4.10 Development of the liquidus temperature at different time instants: (a) t = 10 s, (b) t = 100 s, (c) t = 500 s and (d) t = 1000 s for an initial mass concentration field of 10% and dp = 5 nm. (a) (b) (c) (d) 114 (a) (b) (c) (d) Figure 4.11 Development of solidus temperature at different time instants (a) t = 10 s, (b) t = 100 s, (c) t = 500 s and (d) t = 1000 s for an initial concentration field of 10% and dp = 5 nm. 115 (a) (b) (c) (d) Figure 4.12 Development of the liquidus temperature at different time instants: (a) t = 10 s (b) t = 100 s (c) t = 500 s and (d) t = 1000 s for an initial mass concentration field of 10% and dp = 2 nm. (a) 116 Figure 4.13 Development of solidus temperature at different time instants (a) t = 10 s, (b) t = 100 s, (c) t = 500 s and (d) t = 1000 s for an initial concentration field of 10% and dp = 2 nm. (a) (b) (c) (d) 117 Figure 4.14 Transient development of the nanoparticle concentration profiles at x = 0.5 H and along the y-axis for dp = 5 nm. (a) 118 Figure 4.15 Transient development of the concentration profiles along the y-axis for dp = 2 nm and at (a) x = 0.0043 m (b) x = 0.005 m and (c) x = 0.0055 m. 119 Figure 4.16 Solid-liquid fraction fields at t = 1000 s for: (a) w? = 0, (b) w? = 10%, dp = 5 nm and (c) w? = 10%, dp = 2 nm. (a) (b) (c) 120 Chapter 5 Thermal-Solutal Convection During Solidification of Colloidal Suspensions In this chapter, the model that was developed in Chapter 4 will be used to study thermal- solutal convection during solidification of colloidal suspensions driven by the density difference due to the temperature and concentration variations. The suspensions will be solidified unidirectionally away from the left vertical side of a square cavity. Copper and alumina particles will be used with two different sizes 5 nm and 2 nm, and two cold side temperatures (271 and 268 K) will be utilized. The effect of the rejection rate of the particles, as well as the effect of the mass fraction of the particles, will be investigated and discussed. 5.1 Introduction: Solidification of multi-component mixtures is accompanied with rejection of the solute, which is heavier or lighter than the solvent. This process creates density differences due to the concentration and temperature differences across the mixture, which is the driving force for the creation of thermo-solutal convection currents. The convection currents are responsible for the re-emergence of the solid phase with a different solute composition than the initial melt. Thermo-solutal convection is important in many applications such as metallurgy (Bodeu et al., 2008), oceanography (Turner, 1974), and enhancement of colloidal transport (Selva et al., 2012). An excellent review that explains thoroughly the physics related to the thermo-solutal convection as well as its applications was contributed by Turner (1985). In this brief introduction, we will only focus on the literature on thermo-solutal convection related to the solidification of the 121 binary mixtures. Beckermann and Viskanta (1988) were the first researchers who predicted thermo-solutal transport during a solidification process. They used a hypereutectic solution of water and ammonia, in which the lighter water would be rejected out of the crystallizing phase. They used the one-fluid mixture model as their numerical approach and conducted experimental tests as well. Their results were qualitatively compared with their experiments. They observed that a clockwise rotating flow cell developed in the melt due to thermal convection. However, as time proceeds, a diffusive interface between the top lighter water-rich melt, and the bottom heavier ammonia was developed. A clockwise rotating cell developed as well in the upper layer. Moreover, due to the high water content of the upper layer, the freezing temperature was significantly lowered, which causes the solidification process to be terminated, and the double- diffusive interface vanished, and the melt returns to its initial concentration. Jarvis and Huppert (1995) studied numerically the process of a binary alloy being solidified unidirectionally away from the vertical side. They used two alloys: one that rejects heavier solute than the solvent, and another which rejects lighter solute than the solvent. They found that if the compositionally- heavy solute released from the crystallizing phase, the flow will be driven downward by the thermal and solutal buoyancy currents. However, if the compositionally light fluid is released, there is a possibility of upward-moving convection to emerge. The authors found that the fluid structure that resulted during the solidification process depended on the product of the LeB (Lewis number) and the ratio of the thermal and solutal Rayleigh numbers (?). For the case of ?*LeB1/3>=1, the flow will be thermally-dominated and will exhibit a unidirectional downward flow, whereas for ?*LeB <= 1, the compositional effects overcome thermal effects and unidirectional upper flow will occur. Finally, for the case, Le-1<= ?<= LeB-1/3 , the boundary layer exhibits counterclockwise flow, with the fluid next to the boundary rising, and the fluid further 122 away sinking. The solidification of a binary mixture from the bottom, or the top has been studied extensively experimentally and analytically such in (Wroster, 1986, Kerr et al., 1990, Wettlaufer et al., 1997, Peppin et al., 2007, Peppin et al., 2008). 5.2 Results and Discussion Using the model presented in Chapter 4, results will be presented for a square cavity similar to that of Figure 4.1 that has four solid walls. The suspension is solidified from the left side to the right side of the cavity, and the suspension will consist of water, copper, and alumina nanoparticles. The current direction for the cooling is selected to promote thermal and solutal convection, and investigate different governing parameters such as segregation coefficient, temperature, and particle size on its development. The initial conditions that are implemented are the following: .0tf o r1.0,K274T,0U w i ninin ???? ?? (5.1) The boundary conditions that are imposed are the following: .0,,0,273,0 ,0,0,0,268,0 ,,0,0,0,0 ,0,0,0,0,0 HyHxf or x KTU Hyxf or x KTU HyHxf or yy T U yHxf or yy T U w w w w ???? ? ? ?? ???? ? ? ?? ???? ? ? ? ? ? ? ???? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (5.2) 5.2.1 Development of the liquid fraction and concentration profiles: The transient development of the liquid fraction and the concentration of the particles for the case of copper-water colloidal suspension with a particle size of dp= 2 nm, w? = 10 wt%, k0 123 =0.1, and CT = 268 K are presented in Figures 5.1 and 5.2, respectively. To provide an in-depth understanding of the simultaneous development of the flow field, superimposed liquid fraction field and the instantaneous streamlines at the identified time instants are plotted in Figure 5.1. At the early stage of the solidification process (t=10 s), the solid-liquid interface next to the cold left wall is planar in shape as shown in Figure 5.1a. At the same time, the flow consists of two counter-rotating cells with nearly equal strengths. The counterclockwise rotating cell is observed next to the solid-liquid interface, while the clockwise cell is next to the hot right wall. As time proceeds to t = 100 s (Fig. 5.1b), the appearance of small cell pockets on the interface in the lower part of the interface due to the rejection of the nanoparticles are observed. Particle rejection initiates the constitutional supercooling phenomenon that destabilizes the interface as discussed by El Hasadi and Khodadadi (2013). The strength of the counterclockwise rotating cell and the area that it occupies increased significantly, while the opposite trend was observed for the clock-wise rotating cell with its strength reduced considerably and its coverage area only occupying a small area near the upper-right corner of part the cavity. As time proceeds further as in Figure 5.1c and Figure 5.1d, the interface changes its morphology by the clear appearance of cells in the lower part of the cavity, and the distinct growth of the mushy zone in the lower part of the cavity. At these later times, t= 500 and 1000 s, the liquid melt consists only of a single counterclockwise rotating cell. Here, the crystallized phase rejects particles that are 8.5 times heavier than that of the liquid. In effect, no upward flow convection will emerge. The observations of the current investigation are identical to those of Jarvis and Huppert (1995). The thermo-solutal convection will contribute to the downward flow near the solid-liquid interface, and the flow will be dominated mostly by thermal convection. The transient development of the concentration field for the same working conditions as 124 Figure 5.1 is shown in Figure 5.2. For t = 10 s (Fig 5.2a), a thin zone of the rejected particles is created in front of the solid-liquid interface. However, for t = 100 s (Fig. 5.2b), the zone of the rejected particles has increased in thickness. In addition, some particles have diffused away from the interface. Meanwhile, a zone that is depleted from the particles is formed in the upper left corner of the cavity. However, for t = 500 s, and 1000 s (Figs. 5.2c and 5.2d), the zone of the rejected particles is broken into striations, and the particles are segregated within the channels that are formed between the neighboring dendritic cells of the solid-liquid interface. Some of the particles moved upwards to the hot side of the cavity, due to convection and diffusion, especially because of its high thermal conductivity, which helps the surrounded melt to capture heat, increase its temperature, and become lighter than its surroundings. However, it is intriguing that in the solid phase there are small sites at the bottom of the cavity with high mass fraction of nanoparticles ( w? ? 20 wt%) that indicate that the particles were engulfed by the crystallized phase. One might explain this phenomenon to be due to reduced diffusion of the particles in that specific area so the particles could not escape the sweeping solid-liquid interface. To investigate the effect of the density of the particle on the solidification process of colloidal suspensions, the suspension of water-alumina was investigated under the same conditions as the copper-water suspensions described in Figures 5.1 and 5.2. Alumina particles (?p = 3880 kg/m3) are lighter than copper particles (?p = 8954 kg/m3) by a factor of two and its pertinent properties are listed in Table 5.1. Moreover, alumina is widely used in solutions for preparing freeze casting materials. Figure 5.3 shows a comparison between the superimposed liquid fraction and streamlines for copper and alumina colloidal suspensions at t= 1000 s. The respective solid-liquid interfaces at time t= 1000 s for both cases are almost indistinguishable. However, the dendritic cells were clearly observed and developed for the copper-water 125 suspension (Fig. 5.3a) compared to the alumina suspension (Fig. 5.3b). This resulted because in the copper-water suspension, the solidus and liquidus temperatures are altered more than in the alumina-suspension. For the copper suspensions, the mass fraction of the particles rejected out of the crystalline phase is higher than that for the alumina suspensions. The streamlines for the case of copper-water suspension are packed near the vortex center, whereas for alumina suspension, they were more concentrated near the solid-liquid interface. Higher amounts of particles are accumulated at the lower part of the cavity, at the region between the mushy zone and the crystallizing zone for the case of water-copper suspension (Fig. 5.4a) compared to that for the water-alumina suspension (Fig. 5.4b). This is because copper nanoparticles are heavy, and spread slowly from the regions that are accumulated in front of the solid-liquid interface. Furthermore, the solid phase resulting from the freezing of the alumina nanoparticles contains a higher amount of particles than that for the case of the copper-water suspension. For a better illustration of the effect of the changing of the properties of the nanoparticles on the solidification process, the profiles of the vertical velocity (uy) and the mass fraction of the nanoparticles will be plotted along the horizontal mid-plane (y= H/2) at t = 1000 s. The vertical velocity profiles for alumina and copper suspensions are shown in Fig. 5.5. On the left side of the cavity (x/H<0.3), velocity is zero since this part of the cavity is occupied by the crystallized phase. However, moving away from the fully solidified region, velocity deviates from zero and exhibits small fluctuations while it is a negative value. This region (0.3>x/H>0.4) represents the mushy zone where the liquid coexists with the solid. For the case considered in the current investigation, the intensity of the convention in the mushy region is limited, and that reflects from the small values of the velocity. Away from the mushy area, a sudden increase in the negative value of the vertical velocity (downward flow) is observed and this can be clearly seen 126 by the local minima in Figure 5.5. The copper-water suspension exhibits higher velocity values in the negative direction than the alumina-water suspensions. This can be due to the higher density of copper with respect to that of the alumina, and thus creates a heavier interfacial melt because of the particle's rejection from the solid phase, and as a result, flowing faster downward under the effects of gravity. After the flow attained a local minimum, the velocity direction is shifted from moving downward to moving upwards. This is because the fluid is approaching the hot side of the cavity, getting warmer and lighter until it reaches the local maxim next to the hot wall. The alumina water suspension exhibit higher upward velocities than the copper suspension because it is lighter. The corresponding concentration profiles are shown in Figure 5.6. The concentration is lower than the initial value of 0.1 for most of the crystallized phase for both suspensions considered. However, the solid phase resulted from the freezing of the alumina suspension consists of higher content of nanoparticles than the one that resulted from the freezing of copper suspension even though both suspensions are simulated with the same segregation coefficient. Further investigation is needed in order to elucidate the reasons for this behavior. The concentration attains a local maximum value at the border between the mushy zone and the liquid melt. In the melt, the concentration attains its initial value and remains constant (0.5 0.1 of the segregation coefficient, the interface takes a planar shape, and a thin mushy zone is created in front of the solid-liquid interface. This is due to the low rejection rate of the particles which cannot alter the solidus and liquidus temperatures. In effect, constitutional supercooling is not promoted. However, for the value of k0 = 0.1 (Fig. 5.7d), the interface has transformed from a planar shape to a dendritic one, especially in the part of the interface that is near the lower part of the cavity. In addition, the thickness of the mushy zone is significantly increased due to the increased particles rejected into the liquid melt, which helps to strengthen the constitutional supercooling that is responsible for creating the dendrites. The flow in the liquid melt is characterized by a single counterclockwise rotating cell, which occupies the whole liquid melt. For k0 = 1.0, 0.5, and 0.3 (Figs. 5.7 a, b, and c), the center of the vortex is nearly at the center of the liquid melt. However, for the k0 = 0.1 case, the center of the vortex has moved upward toward the upper-right corner of the cavity. This is due to the increase of the concentration of the particles which increases the intensity of the downward flow and the development of the mushy zone, which helps to push the vortex upward. To further quantify the impact of the segregation coefficient on convection within the melt for the cases presented in Figure 5.7 further, the non-dimensional maximum vorticity is defined as: 0.1max, max, ? ? k kr ??? (5.3) 128 Figure 5.8 shows that the strength of the vortex depends non-monotonically on the segregation coefficient at time of t = 1000 s. The general observation is that particle rejection helps to improve the convection strength if it is compared to the case without particle rejection (k0 = 1.0). 5.2.3 The effect of changing the solidification rate: One of the most crucial factors that influence the solidification process of mixtures is the freezing rate. In the current study, due to the numerical method selected, the interface velocity (i.e. freezing rate) cannot be controlled explicitly. However, increasing the cold side temperature can promote a slower crystallization rate. Therefore, to investigate the effect of the freezing rate, a comparison between the solidification process of CT = 271 K and 268 K will be made with all other conditions kept constant as in section 5.2.1. The superimposed contours of stream functions and liquid fraction fields are shown in Figure 5.9. Clearly the solid content for the case of CT = 268 K (Fig. 5.9a) is higher than for of CT = 271 K (Fig. 5.9b), because the colloid is solidified faster. Furthermore, the mushy zone was extended much further from the solid-liquid interface for the case of higher cold side temperature if it compared to that of the lower one. This is because the particles had a sufficient time to diffuse away from the interface. An interesting feature of the resulting morphology of the solid-liquid interface is that the distance between the dendritic structures is decreased as the cold temperature was decreased (Figure 5.10) similar to experimental observations for colloidal suspensions (Waschkies et al., 2011). 5.2.4 The effect of the particle size on the thermo-solutal convection In the previous chapter, the effect of the particle size on the morphology of the liquid- solid interface was demonstrated. In the current section, the effect of the particle size on the 129 development of the thermo-solutal convection will be presented. Results are presented for the case of w? = 10 wt% water-copper suspension with dp = 5 nm, and 2 nm, and CT = 268 K. The contours of the superimposed stream function and liquid fraction for different times for the case of dp = 5 nm are shown in Figure 5.11. The solid-liquid interface exhibits a planar morphology throughout the freezing process, the mushy zone becomes significantly thin. This is because the liquidus and solidus temperatures are not significantly affected by the concentration of the particles, which prevented the phenomena of constitutional supercooling to be initiated as shown in Chapter 4. The flow field initially consisted of two counter rotating cells, one rotating clockwise near the solid-liquid interface and the other rotating clockwise near the hot wall. However, as time proceeds, the two cells are merged together creating one cell, with its center located at the center of melt region. A brief comparison with the case of dp = 2 nm indicates that the center of the rotating cell is moved upward if it is compared to the case of dp = 5 nm. Moreover, in order to illustrate the effect of the particle size on the thermal-solutal convection, the ratio between the maximum vorticity for the case of dp = 2 nm to that of dp = 5 nm will be plotted at different time instances, and the definition of is given by Eq. (5.4). nmd nmd p p r 5max , 2max , 1 ? ?? ? ?? (5.4) As shown from Figure 5.12, maximum vorticity strength for dp = 5 nm, and 2 nm cases are nearly similar at the initial stages of the freezing process. However, as the solidification process proceeds in time the strength of the vorticity and thus the flow field is higher for the case of dp = 2 nm if it is compared to that of suspension with 5 nm particle size. The stronger convection currents for the case of suspensions with particle size of 2 nm can be attributed to 130 their higher mass diffusion coefficients, compared to that of dp = 5 nm. Higher diffusion coefficient allowed the particles to diffuse much further from the solid-liquid interface, in shorter times and thus helping the mixing of the cold and hot currents of the suspension much more effectively. The concentration fields for the case of dp = 5 nm is shown in Figure 5.13. A layer of rejected particles is formed in front of the solid-liquid interface. As time proceeds (Fig 5.13 b, c, and d), a highly concentrated layer of particles is accumulated at the bottom side of the solid- liquid interface, the crystallized phase resulted consisting of less mass fraction of particles compared to that of the initial state. However, the particles did not diffuse away from the interface as in the case of dp = 2 nm (Figure 5.2) due to their small mass diffusion coefficient. 5.2.4 The effect of the mass fraction of the particles: The mass fraction of the particles plays a significant role in freezing process of colloidal suspensions since all the thermophysical properties are dependent on mass fraction. In this brief section, the effect of the mass fraction of the particles will be investigated for the case of dp = 2 nm, CT = 268 K, and for w? = 0.05, and 0.1 on the morphology of the solid-liquid interface and thermal-solutal convection. The superimposed stream function and liquid fraction for the cases of w? = 0.05 and 0.1 at t = 1000 s are shown in figure 5.14. The solid-liquid interface for the case of w? = 0.05 consists of fewer numbers of dendritic structures if it is compared to that of w? = 0.1. Also, these dendritic structures are smaller in size for the case of w? = 0.05. It can be reasoned that for lower mass fraction of the particles the tendency to change the liquidus, and solidus temperatures are lower, and thus forming the dendrites are less likely. The thickness of the mushy zone for the case of 131 w? = 0.1 is thicker. The stream function field for case w? = 0.05 consisted of a single cell rotating counterclockwise, with the center of the cell nearly at the center of the square cavity. The intensity of the convection cell for the case of w? = 0.1 is higher than that of the case of w? = 0.05. In order to quantify better the development of the convection for two mass fractions, the ratio of the maximum vorticity for the case of w? = 0.1 to that of w? = 0.05 will be plotted for different time instances as shown in Figure 5.15. The value of 2r? is given as the following: 05.0max , 1.0max , 2 ? ?? w wr ? ? ? ?? (5.5) The value of 2r? is increasing with time until it reaches a local maximum as shown in Figure 5.15. This indicates that thermo-solutal convection for the case of w? = 0.1 is more intense than that of w? = 0.05 and thus a better mixing of the suspension occurs. The value of 2r? is nearly equal to 1.5 at t = 500 sec. However, as time proceeds further the value of 2r? is decreasing, which shows that the increase in the value of 2r? with time is not a monotonic function. The concentration contours of w? = 0.1 and 0.05 are shown in Figure 5.16. The main difference between the two concentration fields is that the particle-depleted zones are more developed the case of w? = 0.1 than for of 0.05. 132 5.3 Summary: The following conclusions can be drawn: 1- The rejection of the particles from the solidifying phase into the melt led an increase in the strength of the downward flow. 2- Reducing the value of the segregation coefficient resulted in the increase of thermal- solutal convection intensity. 3- The crystallized phase that resulted from the freezing of the alumina suspension consisted of a higher concentration of particles than that from the solidification of the copper suspension. 4- The crystallized phase that resulted from the lower rate of freezing consisted of a lower concentration of particles than with a higher rate of freezing. 5- As the particle size was reduced, the intensity of thermal-solutal convection as indicated by the variation of the maximum vorticity was increased. 6- As the mass fraction of the particles increased, the intensity of the thermal-solutal convection increased. 133 Table 5.1 Thermophysical and transport properties for the solvent and the nanoparticles Water Alumina nanoparticles Density 997.1 kg/m3 3880 kg/m3 Viscosity 8.9x10-4 Pa s - Specific Heat 4179 J/kg K 729 J/kg K Thermal Conductivity 0.6 W/m K 42.34 W/m K Thermal Expansion Coefficient 2.1x10-4 K-1 5.4x10-6 K-1 Heat of Fusion 3.35x105 J/kg - 134 Figure 5.1 Development of the flow field (shown by the streamlines) superimposed on the contours of the liquid fraction for Ra = 7.44x104, and copper nanoparticles of dp = 2 nm at (a) t = 10 s, (b) t = 100 s, (c) t = 500 s, and (d) t = 1000 s. (a) (b) (c) (d) 135 (a) (b) (c) (d) Figure 5.2 Development of the concentration field for with Ra = 7.44x104, and copper nanoparticles of dp = 2 nm at different time instants (a) t = 10 s, (b) t= 100 s, (c) t = 500 s, and (d) t = 1000 s. 136 Figure 5.3 Superimposed contours of the velocity streamlines and the liquid fraction at t = 1000 s for (a) water-copper (Ra = 7.44x104) and dp = 2 nm and (b) water-alumina suspensions (Ra = 4.36x104) and dp = 2 nm. (a) (b) 137 Figure 5.4 Concentration field at t = 1000s for (a) water-copper, Ra = 7.44x104, and dp = 2 nm and (b) for water-alumina suspensions, Ra = 4.36x104. (a) (b) 138 Figure 2. 7Figure 5.5 Variation of the vertical velocity (uy m/s) profiles at y = 0.5H and along the x-axis at t = 1000 s, and dp = 2 nm. 139 Figure 5.6 Variation of the concentration profiles at y = 0.5H and along the x-axis at t = 1000 s and dp = 2 nm for different particle types. 140 (a) (b) (c) (d) Figure 5.7 Superimposed contours of streamlines and liquid fraction, for Ra = 7.44x104 and dp = 2 nm at t = 1000 s for (a) k0 = 1.0, (b) k0 = 0.5, (c) k0 = 0.3, and (d) k0 = 0.1. 141 Figure 5.8 Variation of the non-dimensional maximum vorticity strength with k0, at t = 1000 s for dp = 2 nm. 142 (a) (b) Figure 5.9 Superimposed contours of streamlines and liquid fraction at t = 1000 s for (a) CT = 268 K, Ra =7.44x104 and for dp = 2 nm, and (b) CT = 271 K, Ra = 2.91x104. 143 (a) (b) Figure 5.10 Concentration field at t = 1000 s, and dp = 2 nm (a) CT = 268 K, Ra = 7.44x104 and (b) CT = 271 K, Ra = 2.91x104. 144 (a) (b) (c) (d) Figure 5.11 Superimposed contours of streamlines and liquid fraction, for dp = 5 nm, Ra = 7.44x104, and CT = 268 K for (a) t= 10 s, (b) t = 100 s, (c) t= 500 s, and (d) t= 1000 s. 145 Figure 5.12 Variation of the non-dimensional maximum vorticity strength with time for different particle sizes 146 (a) (b) (c) (d) Figure 5.13 Contours of the concentration field, for dp = 5 nm, Ra = 7.44x104, and CT = 268 K for (a) t= 10 s, (b) t = 100 s, (c) t= 500 s, and (d) t= 1000 s. 147 (a) (b) Figure 5.14 Superimposed contours of streamlines and liquid fraction at t = 1000 s for CT = 268 K, and dp = 2 nm (a) w? = 0.05, and (b) w? = 0.1, Ra =7.44x104. 148 Figure 5.15 Variation of the non-dimensional maximum vorticity strength with time for different mass fractions of the particles. 149 (a) (b) Figure 5.16 Concentration field contours at t = 1000 s for CT = 268 K, and dp = 2 nm (a) w? = 0.05, and (b) w? = 0.1, Ra =7.44x104. 150 Chapter 6 Conclusions and Perspectives In this chapter, final remarks will be presented to conclude this dissertation. Future work and suggestions will be given for the enhancement of the current prediction methods of the solidification of colloidal suspensions. 6.1 Concluding Remarks: Based on the analytical and numerical results discussed in the previous chapters, the following concluding remarks can be drawn: For the ideal case of no particle rejection out of the solid phase, increasing the thermal diffusivity by adding high conductive particles does not ensure the expedited movement of the solid-liquid interface compared to that of the pure solvent. Since, the equilibrium freezing temperature is reduced as the mass fraction of the particles is increased. At equilibrium, solidification of the colloidal suspensions with particle rejection was controlled by two mechanisms, i.e. thermal and solutal transport. The parameter that controls which mechanism will dominate the other is the transitional segregation coefficient. The transitional segregation coefficient increased with mass fraction of the particles and decreased with the increase in the supercooling and the decrease of the particle size of the suspension. Movement of the interface was decelerated as the value of the segregation coefficient was increased, and the size of the particles was reduced. 151 Reducing the particle size helped to destabilize the solid-liquid interface. For the particle size of 5 nm, the interface has a planar morphology, although for particle size of 2 nm, the interface switched from planar to dendritic morphology because of the constitutional supercooling phenomena for water-copper suspensions. For a planar interface, the rejected particles accumulated in front of the solid-liquid interface, and their concentration increased with time. However, for an unstable interface, the rejected particles accumulated in the space between the dendritic cells, and the dendrite represents the location where the concentration of the particles was the lowest. The rejection of the particles increases the strength of the thermal-solutal convection due to the increase of the downward flow. As the mass fraction of the particles increased, and the value of the segregation decreased, the strength of the thermal-solutal convection was increased. 6.2 Suggestions for Future Work: The following suggestions will expand the current model: I. Use a variable temperature boundary condition at the cold side, to simulate the effect of variable interface velocity on the morphology of the solid-liquid interface and the resulted microstructure. II. Develop a method to characterize the morphology of the solid-liquid interface, and identify the important parameters. 152 III. Use the current model to develop adaptable NePCM that can respond to the external stimuli by an optimum way. IV. Improve the current numerical model to account for the non-equilibrium segregation coefficient, which changes with the velocity of the solid-liquid interface. V. Implement a phase field model (Provatas and Elder, 2010) to simulate solidification of the colloids for the case that the suspension is initially at a supercooling state. In addition, the phase field method had the ability to capture the morphology of the interface more accurately. VI. One of the key elements for a successful simulation of the freezing of colloidal suspension is the accurate prediction of the segregation coefficient. However, at this time, the segregation coefficient is either estimated or calculated from limited thermodynamic relations. However, a promising route for predicting the segregation coefficient is the phase field crystal method that now is used for binary alloy solidification (Humadi et al., 2013). The phase crystal method, which is based on the density functional theory, can numerically simulate the freezing process at an atomic scale and diffusive time. 153 Bibliography ? Acem, Z., Lopez, J., and Palomo Del Barrio, E., 2010, ?KNO3/NaNO3 graphite materials for thermal energy storage at high temperature: Part I. elaboration methods and thermal properties?. Applied Thermal Engineering, 30, pp.1580-1585. ? Azouni, M. A., and Casses, P., 1998, ?Thermophysical Properties Effects on Segregation during Solidification,? Adv. Colloid Interface Sci., 75, pp. 83-106. ? Beckermann, C., and Viskanta, R., 1988, ?Double-Diffusive Convection during Dendritic Solidification of Binary Mixture,? PCH Physicochemical Hydrodynamics, 10, pp. 195- 213. ? Bennon, W. D., and Incropera, F. P., 1987, ?A Continuum Model for Momentum, Heat and Species Transport in Binary Solid-Liquid Phase Change Systems, I Model Formulation,? International Journal of Heat and Mass Transfer, 30, pp. 2161-2170. ? Boden, S., Eckert, S., Willers, B., and Gerbeth, G., 2008, ?X-ray radioscopic visualization of the solutal convection during solidification of a Ga-30%wt Pct in alloy?, Metallurgical and Materials Transactions A, 39, pp. 613-623. ? Buongiorno, J., 2006, ?Convective Transport in Nanofluids,? Journal of Heat Transfer-T ASME, 128, pp. 240-251. 154 ? Catalina, A. V., Mukherjee, A., and Stefanescu, D. M., 2000, ?A Dynamic Model for the Interaction between a Solid Particle and an Advancing Solid/Liquid Interface,? Metallurgical and Materials Transactions A, 31, pp. 2559-2568. ? Chang, A, Dantzig, J., Derr, B. T., and Hubel, A., 2007, ?Modeling the Interaction of Biological Cells with Solidifying Interface,? Journal of Computational Physics, 226, pp. 1808-1829. ? Chiareli, A.O.P., Huppert, H. E., and Worster, M. G., 1999, ?Segregation and flow during the solidification of alloys?, Journal of Crystal Growth, 139, pp. 134-146. ? Cisse, J., and Bolling, G. F., 1971, ?A study of the trapping and rejection of insoluble particles during freezing of water?, Journal of Crystal Growth, 10, pp. 67-76. ? Corte A.E., 1963, ?Particle sorting by repeated freezing and thawing?, Science, 142, 499. ? Corte, A. E., 1962, ?Vertical migration of particles in front of a moving freezing plane?, Journal of Geophysical Research, 67, pp. 1085-1090. ? Crank, J., 1987, Free and moving boundary problems. Oxford University Press, USA. ? Cui, Y., Liu, C., Hu, S., and Yu, X., 2011, ?The experimental exploration of carbon nanofiber and carbon nanotube additives on thermal behavior of phase change materials?, Solar Energy Materials and Solar Cells, 95, Apr., pp. 1208?1212. ? Dantzig, J. A. and Rappaz, M., 2009,?Solidification?, EPFL Press, 620 pages. ? Dash, J. G., Fu, H., and Wettlaufer, J. S, 1995, ?The premelting of ice and its environmental consequences?, Reports on Progress in Physics, 58, pp. 155-167. ? Dash, J. G., Rampel, A. W., and Wettlaufer, J. S., 2006, ? The physics of premelted ice and its geophysical consequences?, Reviews of Modern Physics, 78, pp. 695-741. ? Deville, S., Maire, E., Lasalle, A., Bogner, A., Gauthier, C., Leloup, J., and Gaizard, 155 C., 2009c, ?In situ X-ray radiography and tomography observations of the solidification of aqueous alumina particles suspensions Part II: steady state?, Journal of the American Ceramic Society, 92, pp. 2497-2503. ? Deville, S., Maire, E., Lasalle, A., Bogner, A., Gauthier, C., Leloup, J., and Gaizard, C., 2009b, ?In situ X-ray radiography and tomography observations of the solidification of aqueous alumina particles suspensions Part I: initial instants?, Journal of the American Ceramic Society, 92, pp. 2497-2503. ? Deville, S., Marie, E., Bernad-Granger, G., Lasalle, A., Bogner, A., Gauthier, C., and Guzard, C., 2009a, ?Metastable and unstable cellular solidification of colloidal suspensions?, Nature Materials, 8, pp. 966-972. ? Deville, S., Saiz, E., and Tomsia, A., 2007, ?Ice-templated porous alumina structures?, Acta Materialia, 55, pp. 1965?1974. ? Elliott, J., and Peppin, S., 2011, ?Particle trapping and banding in rapid colloidal solidification?, Physical Review Letters, 107, pp. 1?5. ? Elmer, J. W., Aziz, M. J., Tanner, L. E., Smith, P. M., and Wall, M. A., 1994, ? Formation of bands of ultrafine beryllium particles during rapid solidification of Al-Be alloys: modeling and direct observations. Acta Metallurgica et Materialia, 42, pp. 1065- 1080. ? Fan, L., and Khodadadi, J. M., 2012, ?A theoretical and experimental investigation of unidirectional freezing of nanoparticle-enhanced phase change materials? Journal of Heat Transfer, 133, 092301. 156 ? Garvin, J. W., and Udaykumar, H. S., 2003a, ?Particle-Solidification Front Dynamics Using a Fully Coupled Approach, Part I: Methodology,? Journal of Crystal Growth. 252, pp. 451-466. ? Garvin, J. W., and Udaykumar, H. S., 2003b, ?Particle-Solidification Front Dynamics using a Fully Coupled Approach, Part II: Comparison of Drag Expressions,? Journal of Crystal Growth. 252, pp. 467-474. ? Garvin, J. W., and Udaykumar, H. S., 2005, ?Drag on a Ceramic Particle being Pushed by Metallic Solidification Front,? Journal of Crystal Growth, 276, pp. 275-280. ? Halde, R., 1980, ?Concentration of impurities by progressive freezing?, Water Research, 14, pp. 575-580. ? Ho, C., and Gao, J., 2009, ?Preparation and thermophysical properties of nanoparticle-in- paraffin emulsion as phase change material?, International Communications in Heat and Mass Transfer, 36,pp. 467?470. ? Humadi, H., Hoyt, J. J., and Provatas, N., 2013, ?Phase-filed-crystal study of solute trapping?, Physical Review E, 87, 022404. ? Ishiguro, H., and Rubinsky, B., 1994, ?Mechanical interactions between ice crystals and red blood cells during directional solidification?, Cryobiology, 31, pp. 483-500. ? Jarvise, R. A., and Huppert, H. E., 1995, ?Solidification of a binary alloy of variable viscosity from vertical boundary?, Journal of Fluid Mechanics, 303, pp. 103-132. ? Kerr, R. C., Wood, A. W., Worster, M. G., and Huppert, H. E., 1990, ? Solidification of an alloy cooled from above part 1. Equilibrium growth?, Journal of Fluid Mechanics, 216, pp. 323-342. ? Khodadadi, J. M., and Hosseinizadeh, S. F., 2007, ?Nanoparticle-enhanced phase change 157 materials (NEPCM) with great potential for improved thermal energy storage?, International Communications in Heat and Mass Transfer, 34, pp. 534?543. ? Khodadadi, J. M., Fan, L., and Babaei, H., 2013, ?Thermal conductivity enhancement of nanostructure-based colloidal suspensions utilized as phase change materials for thermal energy storage: A review?, Renewable and Sustainable Reviews, 24, pp. 418- 444. ? Krober, C., and Rau, G., 1985, ?Interaction of Particles and Moving Ice-Liquid Interface,? Journal of Crystal Growth, 72, pp. 649-662. ? Lasalle, A., Guizard, C., Maire, E., Adrien, J., and Deville, S., 2012, ?Particle redistribution and structural defects development during ice templating?, Acta Materialia, 60, pp. 4594-4603. ? Mashl, S. J., Flores, R.A., and Trivedi, R., 1996, ?Dynamics of solidification in 2% corn strach-water mixtures: Effect of variations in freezing rate on product homogeneity?, Journal of Food Science, 61, pp.760-765. ? Nakagawa, K., Thongprachan, N., Charinpanitkul, T., and Tantha-Panichakoou, W., 2010, ?Ice crystal formation in the carbon nanotube suspension: A modeling approach?, Chemical Engineering Science, 65, pp. 1438-1451. ? Peppin, S. S. L., Aussillous, P., Hupport, H. E., and Wroster, M. G., 2007, ?Steady-state mushy layers: experiments and theory?, Journal of Fluid Mechanics, 570, pp. 69-78. ? Peppin, S. S. L., Elliott, J. A. W., and Worster, M. G., 2006, ?Solidification of colloidal suspensions?, Journal of Fluid Mechanics, 554, pp. 147-166. ? Peppin, S. S., Worster, M. G., and Wettlaufer, J. S., 2007, ?Morphological Instability in Freezing Colloidal Suspensions,? Proceedings of Royal Society, 463, pp. 723-733. 158 ? Peppin, S., Wettlaufer, J., and Worster, M., 2008, ?Experimental verification of morphological instability in freezing aqueous colloidal suspensions?, Physical Review Letters, 100(23), June, pp. 1?4. ? Provatas, N., and Elder, K., 2010, ?Phase-field methods for materials science and engineering?, Wiley-Vch. ? Rempel, A., and Worster, M., 1999, ?The interaction between a particle and an advancing solidification front?, Journal of Crystal Growth, 205(3), pp. 427?440. ? Rempel, A., and Worster, M., 2001, ?Particle trapping at an advancing solidification front with interfacial- curvature effects?, Journal of Crystal Growth, 223(3), pp. 420?432. ? Sanusi, O., Warzoha, R., and Fleischer, A. S., 2011, ?Energy storage and solidification of paraffin phase change material embedded with graphite nanofibers?, International Journal of Heat and Mass Transfer, 54, pp. 4429?4436. ? Selva, B., Daubersies, L., and Salmon, J-B., 2012,? Solutal convection in confined geometries: Enhancement of colloidal transport?, Physical Review Letters, 108, 198303. ? Shangguan, S., Ahuja, S., and Stefanescu, D. M., 1989, ?An Analytical Model for the Interaction between an Insoluble Particle and Advancing Solid/Liquid Interface,? Metallurgical and Materials Transactions A, 23, pp. 669-680. ? Shin, D., and Banerjee, D., 2011a, ?Enhancement of specific heat capacity of high- temperature silica-nanofluids synthesized in alkali chloride salt eutectics for solar thermal-energy storage applications?, International Journal of Heat and Mass Transfer, 54, pp. 1064?1070. ? Shin, D., and Banerjee, D., 2011b, ?Enhanced specific heat of silica nanofluid?, Journal of Heat Transfer, 133, 024501. 159 ? Spannuth, M., 2010, ?Structure and dynamics in freezing and frozen colloidal suspensions from direct observation and X-ray scattering?, Ph.D. thesis, Yale University. ? Stan, C., Schneider, G. F., Shevkoplyas, S. S., Hashimoto, M., Ibanescu, M., Wiley, B. J., and Whitesides, G. M., 2009, ? A microfluidic apparatus for the study of ice nucleation in supercooled water drops?, Lab on Chip, 16, pp. 2293-2305. ? Turner, J. S., 1974, ?Double-diffusive phenomena?, Annual Review of Fluid Mechanics, 6, pp. 37-54. ? Turner, J. S., 1985, ?Multicomponent convection?, Annual Review of Fluid Mechanics, 17, pp. 11-44. ? Uhlmann, D. R., Chalmers, B., and Jackson, K. A., 1964, ?Interaction between Particle and Solid-Liquid Interface,? Journal of Applied Physics, 35, pp. 2986-2993. ? Voller, V. R., and Prakash, C., 1987, ?Fixed Grid Numerical Modeling Methodology for Convection-Diffusion Mushy Region Phase Change Problems,? International Journal of Heat and Mass Transfer, 30, pp. 1709-1719. ? Voller, V., 2008, ?A numerical method for the Rubinstein binary-alloy problem in the presence of an under-cooled liquid?, International Journal of Heat and Mass Transfer, 51, pp. 696?706. ? Wakao, N., and Kaguei, S., 1982, Heat and Mass Transfer in Packed Beds, Gordon and Breach Science Publishers, New York, NY, pp. 175-205. ? Wang, N., Zhang, X. R., Zhu, D. S., and Gao, J. W., 2011, ?The investigation of thermal conductivity and energy storage properties of graphite/paraffin composites?, Journal of Thermal Analysis and Calorimetry, 107, Mar., pp. 949?954. 160 ? Waschkics, T., Oberacker, R. and Hoffman, M. J., 2009, ? Control of lamellae spacing during freeze casting of ceramics using double-side cooling as a novel processing route?, Journal of the American Ceramic Society, 92, pp. 79-84. ? Waschkies, T., Oberacker, R., and Hoffmann, M., 2011, ?Investigation of structure formation during freeze-casting from very slow to very fast solidification velocities?. Acta Materialia, 59, pp. 5135?5145. ? Watanabe, K., 2002. ?Relationship between growth rate and supercooling in the formation of ice lenses in a glass powder?. Journal of Crystal Growth, 237-239, pp. 2194?2198. ? Watanabe, K., Muto, Y., and Mizoguchi, M., 2001, ?Water and Solute Distributions near an Ice Lens in a Glass-Powder Medium Saturated with Sodium Chloride?, Crystal Growth and Design, 1, pp. 207-211. ? Wettlaufer, J. S., Wroster, M. G., and Huppert, H. E., 1997, ?Natural convection during solidification of an alloy from above with application to the evolution of sea ice?, Journal of Fluid Mechanics, 344, pp. 291-316. ? Wolfram Mathematica 7 Documentation center, http://reference.wolfram.com/legacy/v7/. ? Wroster, M. G., 1986, ?Solidification of an alloy from a cooled boundary?, Journal of Fluid Mechanics, 167, pp. 481-501. ? Wu, S., Zhu, D., Li, X., Li, H., and Lei, J., 2009, ?Thermal energy storage behavior of Al2O3H2O nanofluids?, Thermochimica Acta, 483, pp. 73?77. ? Yavari, F., Fard, H. R., Pashayi, K., Rafiee, M. A., Zamiri, A., Yu, Z., Ozisik, R., Borca- tasciuc, T., and Koratkar, N., 2011, ?Enhanced thermal conductivity in a nanostructured 161 phase change composite due to low concentration graphene additives?, Journal of Physical Chemistry, 115, pp. 8753?8758. ? Zheng, R., Gao, J., Wang, J., and Chen, G., 2011, ?Reversible temperature regulation of electrical and thermal conductivity using liquid-solid phase transitions?, Nature Communications, 2, 289.