Development of the solid–liquid interface, distribution of the particle concentration field, as well as the development of thermosolutal convection during solidification of colloidal suspensions in a differentially heated cavity are investigated. The numerical model is based on the one-fluid mixture approach combined with the single-domain enthalpy porosity model for phase change, and it is implemented in fluent software package. The linear dependence of the liquidus and solidus temperatures with the concentration of the nanoparticles was assumed. A colloidal suspension consisting of water and copper or alumina nanoparticles were considered. In the current investigation, the nanoparticle size selected was 5 and 2 nm. The suspension was solidified unidirectionally inside a square differentially heated cavity that was cooled from the left side. It was found that the solid–liquid interface changed its morphology from a planar shape to a dendritic one as the solidification process proceeds in time, due to the constitutional supercooling that resulted from the increased concentration of particles at the solid–liquid interface rejected from the crystalline phase. Initially, the flow 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 resulted in crystallized phase with a higher concentration of particles. If it is compared to that of the solid phase resulted from freezing the copper–water colloidal suspension. Decreasing the segregation coefficient destabilizes the solid–liquid interface and increases the intensity of the convection cell with respect to that of the case of no particle rejection. At slow freezing rates, the resulting crystal phase consisted of lower particle content compared to the case of higher freezing rate.
Solidification of multicomponent 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 thermosolutal convection currents. Those convection currents are one of the mechanisms responsible for the solute redistribution in the liquid melt. Thermosolutal convection is important in many applications such as metallurgy [1–3], oceanography , and enhancement of colloidal transport [5–8]. An excellent review that explains thoroughly the physics related to the thermosolutal convection as well as its applications was contributed by Turner . In this brief introduction, we will only focus on the literature on thermosolutal convection related to the solidification of binary mixtures. Beckermann and Viskanta  were among the first researchers who studied thermosolutal 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 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 caused the solidification process to be terminated, the double-diffusive interface vanished, and the melt returned to its initial concentration. Jarvise and Huppert  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 was 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 establishment of upward-moving convection. The authors found that the convection cells 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 , the flow will be thermally dominated and will exhibit a unidirectional downward flow, whereas for , the compositional effects overcome thermal effects and unidirectional upper flow will occur. Finally, for the case , the boundary layer exhibits counterclockwise flow, with the fluid next to the boundary rising, and the fluid further away sinking. The solidification of a binary mixture from the bottom or the top has been studied extensively experimentally and analytically such in Refs. [12–16].
Most of the studies concerned with solidification of colloidal suspensions are based on experiments due to the lack of theoretical progress on the subject. Solidification of colloidal suspensions have a wide range of applications such as in food-processing , purification of water , preservation of blood cells , creating new composite materials [20,21], and energy storage . Most of the investigations were conducted with a cold bottom side [23–25]. It was observed that the morphology of the interface is changed from planar to cellular to lamellar, and finally to dendritic. They found that the cooling rate, the volume fraction of the particles, and particles size are the key factors that affect the morphology of the resulting microstructure.
However, the theoretical investigations of colloidal suspension solidification are limited to analysis of 1D Stefan type models as in Refs. [26–29]. They found that as the volume fraction of the particles increases and the size of the particles decreases, the solid–liquid interface changes from a planar stable one to an unstable morphology. Recently, El Hasadi and Khodadadi  investigated the solidification process of water–copper colloidal suspension confined in a two-dimensional rectangular cavity, for a single mass fraction of particles of 10% and two particle diameters of 2 nm and 5 nm. They used the one-fluid mixture model, and the suspension was frozen from the bottom. They observed that, for the case of 2 nm, the interface in the late stage of the solidification process developed a dendritic morphology, due to the development of the constitutional supercooling phenomena resulting from rejection of particles in the interdendritic region. Furthermore, convection currents in the interdendritic region were observed.
To the best knowledge of the authors, there is no study in the literature which includes the impact of the segregation coefficient, different particle properties, and the freezing rate on the solidification process of the colloidal suspensions. The effect of these parameters and of the particle size on thermosolutal convection resulting from the freezing process will be investigated as well. In this computational study, representative copper–water and alumina–water suspensions are solidified inside a square cavity, frozen from the vertical side. Different segregation coefficients and various cold side temperatures are investigated.
Since colloidal suspensions can be considered as binary mixtures, their solidification process is accompanied with 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 it. The model that will be used for this study is based on the one-fluid mixture model that has been used extensively for modeling binary alloy solidification such as in Refs. [10,31] which treats the mushy zone as a porous medium. The assumptions used in the current model are as follows:
The solvent is assumed to be a homogeneous material in which the solid and the liquid phases have the same thermophysical properties. The solvent properties are taken to be those of liquid water in both liquid and solid phases, and the density inversion phenomenon was ignored, as a simplification to reduce the nonlinearity of the current system of equations.
Local thermodynamic equilibrium is assumed throughout the simulation.
The solid phase is always stationary.
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.
The thermophysical and transport properties of the dilute () colloid suspension are not constant and change with the concentration of nanoparticles.
Consistently with Assumption 2 of equal densities between phases, convection due solely to change of phase is ignored.
The geometry of the physical model considered in the present investigation is a square cavity (side H = 10 mm) containing the colloidal suspension as shown in Fig. 1. The freezing of the suspension is initiated by lowering the temperature of the left side of the cavity TC below the liquidus temperature corresponding to the initial concentration, while the right vertical side was kept at a constant temperature (TH = 273 K). The remaining two horizontal sides are kept thermally insulated. The initial and boundary conditions implemented for all the cases considered in the current paper are as follows:
Mixture Relations and Effective Thermophysical/Transport Properties.
The values of the segregation coefficient (k0) ranged from k0 = 1.0 (no particle rejection) to 0.1. The value of the segregation coefficient is kept constant throughout each simulation.
The governing equations used in the current investigation have been discretized using the finite volume method. A uniform staggered grid was used, and the control volumes had equal sides in the y and x directions. The velocity and pressure fields are coupled using the semi-implicit method for pressure-linked equations (SIMPLE) . The second-order upwind scheme  was used to handle the coupling of convection and diffusion terms in the momentum, energy, and species equations. The phase-change process of solidification was simulated using the enthalpy porosity numerical method. The convergence criterion used required that the residuals after the end of each time step for the momentum equations be less than 10−5, whereas for energy and species equations, the value was less than 10−7. The time step was set to 0.001 for the first 1000 time steps, then increased to 0.05 for the next 1000 steps, and finally, it was assigned the value of 0.1 for the rest of the simulation time. The commercial fluid solver fluent has been used to implement numerically the model discussed above. The sets of the equations are interlinked through the source term or through the thermophysical and transport properties for which special user defined functions were developed to handle their variation with the volume fraction of the particles. A grid independence study was conducted using three grid configurations 50 × 50, 100 × 100, and 150 × 150; the results of the grid independency are shown in Fig. 2. The configurations with 100 × 100 grids, have been selected for accuracy and reduced time. For the purpose of benchmarking the adopted model, the results have been compared with those of Hannoun et al. . The comparison was satisfactory and the details of the current model are presented in Ref. .
Results and Discussion
In this section, the results of the numerical investigations will be discussed. Emphasis will be given to the effect of the segregation coefficient, the freezing rate, and particle size. The properties of the nanoparticles and the solvent are listed in Table 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, = 10 wt.%, k0 = 0.1, and = 268 K are presented in Figs. 4 and 5, 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 Fig. 3. 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 Fig. 3(a). 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. 3(b)), the appearance of small cell pockets on the interface in the lower part of the interface, due to the rejection of the nanoparticles, is observed. Particle rejection initiates the constitutional supercooling phenomenon that destabilizes the interface as discussed by El Hasadi and Khodadadi . Simultaneously, the strength of the counterclockwise rotating cell and the area that it occupies increased significantly, while the opposite trend was observed for the clockwise rotating cell with its strength reducing considerably, and its coverage area only occupying a small area near the upper-right corner of the cavity. As time proceeds further as, in Figs. 3(c) and 3(d), 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 the liquid. These observations match those of Jarvise and Huppert . The thermosolutal convection contributes to the downward flow near the solid–liquid interface, and the flow is dominated mostly by thermal convection.
The transient development of the concentration field for the same working conditions as Fig. 3 is shown in Fig. 4. At t = 10 s (Fig 4(a)), a thin zone of the rejected particles is created in front of the solid–liquid interface. However, at t = 100 s (Fig. 4(b)), the zone of the rejected particles has increased in thickness. In addition, some particles have diffused away from the interface. Meanwhile, a zone depleted from particles is formed in the upper left corner of the cavity. However, at t = 500 s and 1000 s (Figs. 4(c) and 4(d)), 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 upward to the hot side of the cavity, due to convection and diffusion, especially because of its high thermal conductivity, which helps the surrounding 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 ( 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 particle density 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 Figs. 3 and 4. Alumina particles (ρp = 3880 kg/m3) are lighter than copper particles (ρp = 8954 kg/m3) by a factor of 2 and its pertinent properties are listed in Table 1. Alumina is widely used in solutions for preparing freeze casting materials. Figure 5 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 suspension (Fig. 5(a)) compared to the alumina suspension (Fig. 5(b)). 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. More 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. 6(a)) compared to that for the water–alumina suspension (Fig. 6(b)). 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 more particles than for the case of the copper–water suspension.
For a better illustration of the effect different properties of the nanoparticles on the solidification process, the profiles of vertical velocity (uy) and of mass fraction of the nanoparticles is plotted along the horizontal midplane (y = H/2) at t = 1000 s. The vertical velocity profiles for alumina and copper suspensions are shown in Fig. 7. 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 its negative value indicates downward motion of the colloid. 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 convection in the mushy region is limited, and that is reflected in small values of velocity. Away from the mushy area, a sudden increase in downward flow is observed and this can be clearly seen by the local minima in Figure 6. The copper–water suspension exhibits higher downward velocity than the alumina–water suspensions. This can be due to the higher density of copper relative to the alumina, which 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 upward. This is because the fluid is approaching the hot side of the cavity, getting warmer and lighter until it reaches the local maximum next to the hot wall. The alumina water suspension exhibit higher upward velocities than the copper suspension because it is lighter. The corresponding instantaneous concentration profiles are shown in Fig. 8. The local concentration is lower than the initial value of 0.1 for most of the crystallized phase for both suspensions considered. However, the solid phase resulting from freezing of the alumina suspension consists of higher content of nanoparticles than from 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 profile attains a local maximum value at the border between the mushy zone and the liquid melt. Within the melt, the concentration attains its initial value of 0.1 and remains constant (0.5 < x/H < 0.9) until it reaches the region near the hot wall where it increases again due to convection and diffusion.
Effect of Changing the Segregation Coefficient (k0).
To investigate the effect of the rejection rate of the particles on the solidification process, different values of the segregation coefficient were selected ranging from no rejection case (k0 = 1.0) to a high rejection (k0 = 0.1). For no particle rejection, particles that are initially distributed homogeneously in the cavity remain static during the freezing process, similar to Khodadadi and Hosseinizadeh . The suspension selected was a copper–water colloid with a particle size of dp = 2 nm, = 268 K, and = 10%. The superimposed contours of the stream function and the liquid fraction field are shown in Fig. 9 for different values of k0 at t = 1000 s. For values k0 > 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 does not alter the solidus and liquidus temperatures. In effect, constitutional supercooling is not promoted. However, for the value of k0 = 0.1 (Fig. 9(d)), 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 increase particle rejection 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 entire liquid melt. For k0 = 1.0, 0.5, and 0.3 (Figs. 9(a)–9(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 particle concentration, which increases the intensity of the downward flow and the development of the mushy zone, helping to push the vortex upward.
Figure 11 shows that the strength of the vortex depends nonmonotonically on the segregation coefficient at time instant t = 1000 s. The general observation is that particle rejection helps to improve the convection strength as compared to the case without particle rejection (k0 = 1.0). The nonmonotonic behavior observed in Fig. 10 can be explained that by decreasing the value of the segregation coefficient the amount the heavy particles rejected into the liquid is increased which results in heavier liquid, and thus the overall convection strength is deceased. The existence of the maximum point probably is the result of a balance between the thermal and solutal effects on the overall convection.
The Effect of Changing the Solidification Rate.
One of the most crucial factors that influences the solidification process of mixtures is the freezing rate. In the current study, the interface velocity (i.e., freezing rate) cannot be controlled explicitly. However, increasing the cold side temperature can promote a slower crystallization rate. Thus, to investigate the effect of the freezing rate, a comparison between the solidification process of = 271 K and 268 K cases will be made with all other conditions kept same as in Sec. 4.1. The superimposed contours of stream function and liquid fraction fields are shown in Fig. 12. Clearly, the solid content for the case of = 268 K (Fig. 11(a)) is greater than the = 271 K case (Fig. 11(b)), because the colloid is solidified faster. Furthermore, the mushy zone was extended much further from the solid–liquid interface in the case of higher cold side temperature. This is because the particles had sufficient time to diffuse away from the interface. A feature of the solid–liquid interface is that the distance between the dendritic structures is decreased as the cold temperature was decreased (the space between the channels for Fig. 12(a) is less than Fig. 12(b)) similar to experimental observations for colloidal suspensions .
The Effect of the Particle Size on the Thermosolutal Convection.
In the current section, the effect of the particle size on the development of the thermosolutal convection will be presented. Results are presented for the case of = 10 wt.% water–copper suspension with dp = 5 nm and 2 nm, = 268 K, and k0 = 0.1.
The contours of the superimposed stream function and liquid fraction at different times for the case of dp = 5 nm are shown in Fig. 13. The solid–liquid interface exhibits a planar morphology throughout the freezing process and the mushy zone remains 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. 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 merged together creating one cell, with its center located at the center of melt region. For the case of dp = 2 nm, the center of the rotating cell is moved upward 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 of 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 is given by Eq. (18).
As shown in Fig. 14, maximum vorticity strength for dp = 5 nm and 2 nm cases are nearly similar at the initial stages of the freezing process. However, as solidification proceeds in time, the strength of the vorticity and thus the flow field is higher for the case of dp = 2 nm. The stronger convection currents for the of dp = 2 nm, case can be attributed to their higher mass diffusion coefficient. 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 Fig. 15. A layer of rejected particles is formed in front of the solid–liquid interface. As time proceeds (Figs. 17(b)–17(d)), a highly concentrated layer of particles is accumulated at the bottom side of the solid–liquid interface and the resulting crystallized phase has lower mass fraction of particles than initially. However, the particles did not diffuse away from the interface as in the case of dp = 2 nm (Fig. 4) due to their small mass diffusion coefficient.
The Effect of the Local Mass Fraction of Particles.
The mass fraction of particles plays a significant role on freezing process of colloidal suspensions, since all the thermophysical properties depend on mass fraction. In this section, the effect of mass fraction of particles will be investigated for the case of dp = 2 nm, = 268 K, and for = 0.05 and 0.1 and k0 = 0.1.
The superimposed stream function and liquid fraction for the cases of = 0.05 and 0.1 at t = 1000 s are shown in Fig. 16. The solid–liquid interface for the case of = 0.05 consists of fewer dendritic structures than for = 0.1. Also, these dendritic structures are smaller in size for the case of = 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 dendrites is less likely. The thickness of the mushy zone for the case of = 0.1 is thicker than that of = 0.05.
The value of is increasing with time until it reaches a local maximum as shown in Fig. 17. This indicates that thermosolutal convection for the case of = 0.1 is more intense than that of = 0.05 and thus a better mixing of the suspension occurs. The value of is nearly equal to 1.5 at t = 500 s. However, as time proceeds further the value of is decreasing, which shows that the increase in the value of with time is not a monotonic function, because as time proceeds the suspension with = 0.1 becomes heavier than that with 0.05 due to the rejection of the particles.
The concentration contours of = 0.1 and 0.05 are shown in Fig. 18. The main difference between the two concentration fields is that the particle-depleted zones are more developed for the case of = 0.1 than for the case of = 0.05.
In the current paper, the thermophysical properties of liquid water and ice are assumed to be equal in order to simplify the numerical model used. However, the change in the properties will affect the resulted solid–liquid interface, the concentration of the particles, and the resulted thermosolutal convection. Thus, further investigation is needed to illustrate the effect of thermophysical properties of the solid and liquid phases of the solvent on solicitation process of colloidal suspensions
A numerical investigation of the solidification process of colloidal suspensions has been conducted, with the suspension treated as binary mixture. Emphasis was given to the effect of various parameters such as the segregation coefficient, cooling rate, the volume fraction, and particle size of the suspensions on the thermosolutal convection that resulted due to the freezing process. The following conclusions were drawn:
The rejection of particles from the solidifying phase into the melt led to an increase in the strength of the downward flow.
Reducing the value of the segregation coefficient resulted in increased of thermal–solutal convection intensity.
The crystallized phase that resulted from the freezing of the alumina suspension had a higher concentration of particles than that from the solidification of the copper suspension.
The crystallized phase that resulted from lower rate of freezing had a lower concentration of particles than with a higher rate of freezing.
As the particle size was reduced, the intensity of thermal–solutal convection, as indicated by the variation of the maximum vorticity, was increased.
As the mass fraction of the particles increased, the intensity of the thermal–solutal convection increased.
This material is based upon work supported by the US Department of Energy under Award No. DE-SC0002470. The first author acknowledges the support of the Department of Mechanical Engineering at Auburn University through providing a teaching assistantship.
- cp =
specific heat (J/kg K)
- Cm =
porosity constant (kg/m3 s)
- dp =
particle diameter (nm)
- DB =
Brownian diffusion (m2/s)
- DT =
thermophoretic diffusion (m2/s)
unit vector in the y-direction
- H =
length of the cavity (m)
- k =
thermal conductivity (W/m K)
- kB =
Boltzmann constant (J/K)
- k0 =
- L =
latent heat (J/kg)
- ml =
liquidus slope (K)
- Ra =
- p =
- t =
- T =
- Tref =
reference temperature (K)
velocity vector (m/s)
- vp =
volume of the particle (m3)
- x, y =
Cartesian coordinates (m)