Accepted Manuscript Title: A lattice Boltzmann model for solidification of water droplet on cold flat plate Author: Jinjuan Sun, Jianying Gong, Guojun Li PII: DOI: Reference:
S0140-7007(15)00193-0 http://dx.doi.org/doi:10.1016/j.ijrefrig.2015.07.003 JIJR 3083
To appear in:
International Journal of Refrigeration
Received date: Revised date: Accepted date:
15-7-2014 6-6-2015 5-7-2015
Please cite this article as: Jinjuan Sun, Jianying Gong, Guojun Li, A lattice Boltzmann model for solidification of water droplet on cold flat plate, International Journal of Refrigeration (2015), http://dx.doi.org/doi:10.1016/j.ijrefrig.2015.07.003. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A lattice Boltzmann model for solidification of water droplet on cold flat plate Jinjuan Sun, Jianying Gong*, Guojun Li MOE Key Laboratory of Thermo-Fluid Science and Engineering, School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China *Corresponding author, Tel: +86 29 82668728; Tel: +86 29 82668789 Email address:
[email protected]
Highlights 1. A 2-D lattice Boltzmann model for solidification of water droplet is proposed. 2. The influence of cold flat plate temperatures on solidification process is studied. 3. Temperature distributions and solid phase variations in droplet are obtained. 4. The shape of solid phase in freezing can also be predicted.
Page 1 of 21
ABSTRACT: Since the solidification of water droplet is the initial and essential process in the whole process of frosting, a model is developed by the lattice Boltzmann method (LBM) that applies the velocity and temperature distribution functions to investigate the solidification process of water droplet on cold flat plate. The thermal transport and liquid-solid phase transition in present model are both based on the pseudo-potential model combined with the enthalpy formation. By this LB model, the solidification process is simulated in form of temperature and solid phase variations in water droplet on cold flat plate, and the shape of solid phase in freezing can also be predicted. In addition, we apply the present LB model to preliminarily study the frost formation process. Numerical results agree well with our experimental data.
Keywords: Lattice Boltzmann method (LBM); Solidification; Water droplet; Cold flat plate. Nomenclature
u
velocity (m s-1)
cp
heat capacity (J K-1)
V
volume of water droplet (μl)
cs
isothermal sound velocity (m s-1)
v
velocity of air (m s-1)
ei
discrete velocity
x
position
F Fg Ft
total force (N) gravity force (N) fluid-solid interaction force (N)
Greek symbols thermal diffusivity (m2 s-1) α
Fσ
fluid-fluid interaction force (N)
δh
fi
velocity distribution function
δx
latent heat of phase transition (J kg-1) lattice spacing
fs
volume fraction of solid phase
δt
time step
G
fluid-fluid interaction coefficient
θ
contact angle (rad)
Gt
fluid-solid interaction coefficient
υ
kinematic viscosity (m2 s-1)
G(x,x')
Green's function total enthalpy (J) enthalpy value at the end of freezing (J) enthalpy value at the start of freezing (J) Boltzmann constant (J K-1) universal gas constant (J mol-1 K-1)
ρ σ τ, τc Φt φ(x)
density of fluid phase (kg m-3) interfacial energy (kJ m-2) dimensionless relaxation time
ωi
weight coefficient
H He Hs k R
source term function of local densities
RH r
relative humidity droplet radius (mm)
Subscripts and superscripts
s(x+ei)
indicator function of solid phase
a
air
T
temperature (°C)
w
cold flat plate
Page 2 of 21
Ti
temperature distribution function
i
direction in a lattice
t
Time (s)
eq
equilibrium
1. Introduction The frost formation has important effects on many fields, such as refrigeration industry, air-conditioning, aviation, etc. The process of frost formation combined with transient heat and mass transfer is extremely complex, which will happen only if the temperature of cold flat plate is less than zero degrees Celsius and the ambient dew point temperature. The presence of frost will deteriorate refrigeration capacity due to lower heat transfer coefficient and higher pressure drop, and defrosting will lead to higher energy consumption and operation cost. Hayashi et al. (1977) photographed the whole process of frost formation, which was divided largely into three stages: crystal growth period, frost layer growth period and frost layer full growth period. Over the last several decades, research has focused on frost layer growth and frost layer full growth period by physical-mathematical and experimental methods (Lee et al., 1997; Yun et al., 2002; Yang et al., 2006; Matsumoto et al., 2014). As the inception of crystal growth period, the solidification of water droplet has an essential role in the whole process of frosting and attracted more and more attention recently. The study of solidification has been conducted in many fields, such as metal processing (Mortensen and Jin, 1992), solidification of castings (Cantor and O'Reilly, 2002) and energy storage system (Esen and Ayhan, 1996; Esen et al., 1998; Esen, 2000), but few aiming at water droplet solidification in frost formation. Hoke et al. (2000) used a high speed camera to observe the freezing process of droplets on a variety of plates and found that the structure and form of the ice immediately after freezing is plate dependent. Tabakova and Feuillebois (2004) constructed a numerical model by the finite-difference fully implicit scheme which can be used to simulate the solidification and subsequent cooling of a supercooled liquid droplet. Wu et al. (2007) investigated the freezing of supercooled condensate droplets in the frost formation, and the freezing onset time and diameter of condensate droplets were characterized. Bahadur et al. (2011) developed a physics-based model to predict ice formation on cooled superhydrophobic plates resulting from the impact of supercooled water droplets, and this modeling approach analyzed the multiple phenomena influencing
Page 3 of 21
ice formation on superhydrophobic plates. Despite several studies have been conducted on solidification of water droplet on cold plate, both the temperature distribution and phase transition process in droplet are still not entirely clear from experimental and analytical results. Thus, an elaborate numerical simulation considering thermal transport is required to predict the phase transition process and temperature field of such freezing water droplet. The lattice Boltzmann method (LBM) has become a particularly successful scheme to simulate fluid flow and complex physical phenomena (Shan and Chen, 1993; Chen and Doolen, 1998; Yu et al., 2003; Aidun and Clausen, 2010; Sheikholeslami et al., 2014). LBM can also effectively simulate the heat and mass transfer in porous media and phase transition. Jiaung et al. (2001) acted as pioneer in applying the enthalpy equation in LB model to simulate liquid-solid phase transition, and numerical results obtained from the improved LB model agreed well with previous analytical or numerical results in the literatures. Guo et al. (2002) first proposed a LB model considering porosity in the equilibrium distribution function, which was improved by Guo and Zhao (2005) to solve convection heat transfer problems and simulate temperature field in porous media. Semma et al. (2008) adopted LBM for solving melting and solidification problems, which was first validated for natural convection with or without coupling of phase transition. Eshraghi and Felicelli (2012) developed a LB model to solve the heat conduction problem with phase transition, and good accuracy was demonstrated compared to finite element method. LBM has the advantages of simple implementation, capability for simulating highly complex geometries and boundaries, computational efficiency and inherent parallel-processing structure in comparison with traditional numerical methods (Eshraghi and Felicelli, 2012). In addition, owing to its statistical roots, LBM might be a more elaborate method to cover the kinetic aspects of the liquid/solid phase transition, especially the interaction of flow and growth kinetics (Miller, 2001). In present paper, our purpose is to simulate the solidification process of water droplet combined with both phase transition and complex boundaries. However, different LB models have different limitations in scopes and conditions of application (Guo and Zheng, 2009; He et al., 2009). The phase transition model and pseudo-potential model are applied in present paper. The phase transition model based on the enthalpy is unable to solve the influence of the phase transition on the flow. The pseudo-potential model cannot simulate the flow with high viscosity or density ratio. To avoid the
Page 4 of 21
limitation of phase transition model based on the enthalpy, the percolation theory in porous media is applied to analyze the flowing problem in phase transition process. And the pseudo-potential model is suitable in present simulation due to non-high viscosity or density ratio. Therefore, LBM is chosen as a suitable and effective method to simulate the solidification of water droplet. The purpose of this paper is to establish a lattice Boltzmann model to numerically study the solidification process of water droplet on cold flat plate, prove this model is efficient by our experimental results and lay the foundation for investigating the extremely complex process of frost formation. The numerical model and theory of lattice Boltzmann method is discussed in Section 2. Section 3 describes an experimental system for solidification of droplet. The variations of temperature distribution and solid phase in water droplet are obtained and the present LB model is verified in Section 4. Section 5 preliminarily studies the frost formation process, followed by a brief conclusion in Section 6.
2. Numerical model and theory
2.1 Lattice Boltzmann equation for the velocity field Several models based on lattice Boltzmann equation can be used to simulate the multiphase fluid flow. One of the most well-known models is the pseudo-potential model (Shan and Chen, 1994), which has advantages of tracking and separating the interface automatically, high computational efficiency, easiness to handle the complex multiphase flow and so on. Therefore, the present paper chose the pseudo-potential model as the basic method for resolving the velocity field. The lattice Boltzmann equation (LBE) is expressed in the following form: f i x e i t , t t f i x , t
1
[ fi x, t fi
( eq )
x , t ]+ tFi
(1)
where fi is the particle distribution function with discrete velocities ei at position x and time t. We adopted the Bhatnagar-Gross-Krook (BGK) (Guo and Zheng, 2009; He et al., 2009) approximation with two-dimension and nine-velocity LB model. Schematic illustration of the lattice Boltzamnn structure is shown in the Fig. 1. The discrete velocities ei in Eq. (1) are therefore written as:
Page 5 of 21
0 x ei (cos[( i - 1) ], sin[( i - 1) ]) 2 2 t x (cos[(2 i - 1) ], sin[(2 i - 1) ]) 2 t 4 4
i 0, i 1, 2, 3, 4,
(2)
i 5, 6, 7, 8.
where δx is the lattice spacing, δt is the time step, and τ is the dimensionless relaxation time connected with the kinematic viscosity c s ( 0.5) t . 2
The equilibrium distribution function (EDF) fi(eq) in Eq. (1) is described as follows: fi
i [1
( eq )
ei u 2
(e i u )
cs
4
2 cs
2
u
2 2
]
2 cs
(3)
where ρ and u are the density and velocity, respectively, the isothermal sound velocity c s x RT / t with RT 1 / 3 , and the weight coefficients for different directions
ωi is given as follows: i0
4 / 9 i 1 / 9 1 / 3 6
i 1, 2, 3, 4
(4)
i 5, 6, 7 , 8
In the pseudo-potential model, the local fluid density and flow velocity are calculated according to the following equations: 8
(5)
fi
i0
u
1
8
( ei fi F )
(6)
i0
where the total force F Fσ Ft Fg , it contains three forces including the interaction force Fσ between fluids, the interaction force Ft between fluid and solid, and the gravity force Fg. The interaction force Fσ between fluids is written as follows: Fσ ( x ) G ( x , x ') ( x ')( x ' x )
(7)
x'
where the function of local densities ( ) 0 [1 exp( / 0 )] (Shan and Chen,
Page 6 of 21
1994). The Green's function G(x,x') in Eq. (7) is written as follows: | x x ' | x
G G ( x , x ') G / 4 0
2 x
| x x ' |
(8)
otherw ise
where G is the interactive strength coefficient between fluids. The interaction force Ft between fluid and solid (Martys and Chen, 1996) is expressed as follows: Ft ( x ) G t ( x , x ')s ( x e i )( x ' x )
(9)
x'
where s(x+ei) is the indicator function of the solid phase. For the fluid phase s(x ei ) 0
and for the solid phase s ( x e i ) 1 . The function Gt(x,x') in Eq. (9) is
written as follows: G t G t ( x , x ') G t / 4 0
| x x ' | x | x x ' |
2 x
(10)
otherw ise
where Gt is the interactive strength coefficient between fluid and solid. Performing a Chapman-Enskog expansion in the pseudo-potential model, the continuity and momentum equations at the Navier-Stokes level are derived from the LBE. t u t
u 0
uu
(11) p
F
u
(12)
2
2.2. Lattice Boltzmann equation for the temperature field In LB model, the evolution equation of the temperature field based on the phase transition is expressed as follows: Ti x e i t , t t Ti x , t
1
c
[T i x , t T i
( eq )
x , t ] t t
(13)
where Ti is the temperature distribution function and τc is the dimensionless relaxation time. In the D2Q9 lattice, the thermal diffusivity c s2 ( c 0.5) t . The equilibrium temperature distribution function Ti(eq) in Eq. (13) is written as follows:
Page 7 of 21
Ti
( eq )
i T [1
ei u 2
(e i u )
cs
2
4
2 cs
u
2 2
]
(14)
2 cs
where the weight coefficients ωi and the isothermal sound velocity cs are the same as those for the velocity field in Eq. (3), and T is the macroscopic temperature defined as follows: 8
T
T
(15)
i
i0
In the Eq. (13), the discretized form of the source term is established as the following equation: h fs fs '
t i
cp
t
(16)
where δh is the latent heat of phase transition, cp is the heat capacity at constant pressure and fs is the volume fraction of solid phase calculated using the following form: 0 H H fs s Hs He 1
H Hs He H Hs
(17)
H He
where Hs and He represent the enthalpy values at the start and end of the freezing, respectively, and H is the total enthalpy calculated as follows: H c pT f s h
(18)
On the basis of the Chapman-Enskog expansion, we can derive the temperature macroscopic equation at the Navier-Stokes level (Jiaung et al., 2001) from the LBE.
2.3. The evolution of LBM The evolution of LBM includes two independent particle movements of the collision and migration. These two movements alternate with time and the sequence is equivalent, but there are differences for simulation in phase transition region. When external temperature is not considered, particles maintain the dynamic equilibrium. The collision and migration in this state will occur simultaneously. However, once the external temperature is involved in the evolution of phase transition, the particle in the boundary will migrate and transfer temperature into system. Then due to the
Page 8 of 21
temperature change the particle collision occurs. Thus the way of migration before collision is more reasonable. Migration: f i x + e i t , t + t = f i x , t
(19)
Ti x + e i t , t + t = Ti x , t
(20)
'
'
Collision: f i x + e i t , t + t f i x + e i t , t + t '
fi
Ti x + e i t , t + t Ti x + e i t , t + t '
Ti
[ f i x + e i t , t + t '
(21)
x + e i t , t + t ]
( eq )
( eq )
1
1
c
[Ti x + e i t , t + t '
(22)
x + e i t , t + t ] t t
To solve the flowing problem in solid-liquid mixture and solid region, we adopted the percolation theory in porous media (Semma et al., 2008). The migration is calculated by modifying Eq. (19) with the volume fraction of solid phase: f i x + e i t , t + t = (1 f s ) f i x , t f s f i x + e i t , t '
(23)
If f s 0 , the node is totally fluid and the migration of particle is as the same as that in Eq. (19). If f s 1 , the node is totally solid and particle remains stationary. If 0 f s 1 , the node is solid-liquid mixture and only a part of f i x , t is migrated.
2.4. Boundary conditions Appropriate boundary conditions must be implemented through the distribution functions in the LBM, which will significantly influence the stability and accuracy of simulation. There are many different types of boundary conditions such as flat boundary condition, surface boundary condition, pressure boundary condition (Guo and Zheng, 2009) and so on. In LB model, the bounce back condition is applied on the nodes of cold flat plate and solid crystal surfaces, which conforms to the quasi-static hypothesis. The periodic boundary condition is applied on the nodes of air boundaries. The temperature of cold flat plate is set by measured value, and the temperature of humid air is kept constant. For the thermal boundary condition, we employed the
Page 9 of 21
non-equilibrium extrapolation scheme (Guo et al., 2002), which is described as follows: Ti x b , t T i
( eq )
x b , t + Ti ( neq ) x f , t
(24)
where x b is a lattice node on the boundary and x f is its nearest neighboring node. The non-equilibrium state is expressed as follows: Ti
( neq )
x f , t Ti x f , t Ti ( eq ) x f , t
(25)
The temperature distribution of each node is calculated by substituting Eq. (25) into Eq. (24): Ti x b , t Ti
( eq )
x b , t + Ti x f , t Ti ( eq ) x f , t
(26)
The non-equilibrium extrapolation scheme of boundary treatment has higher simulative stability, less computation and wider range of application, which also has two-order accuracy in both time and space.
2.5. Procedure of LB model In the physical model, the liquid droplet is cooled to a temperature below the equilibrium freezing point. The supercooling drives rapid kinetic crystal growth from the crystal nuclei. And an abrupt temperature rise occurs in droplet because of released latent heat from the nucleation. This stage is terminated when the supercooling is exhausted and the droplet has reached its equilibrium freezing temperature. Then the growth of the solid phase is governed by heat transfer rate from the droplet to the point where the droplet liquid is completely frozen. Finally, the solid droplet temperature reduces to a steady state value near that of the cold flat plate. The solidification process includes four stages: pre-cooling or supercooling, recalescence, freezing, and solid cooling or tempering. (Hindmarsh et al., 2003) In present paper, the first stage happens when t < 0 and the last two stages happen when t > 0. The recalescence stage is instantaneous and not considered in present simulation. The initial time of freezing (t = 0) begins at the end of the recalescence stage. In order to model the solidification process with different properties, a FORTRAN90 code is written to solve Lattice Boltzmann equation. Before running procedure, physical variables (such as density, humidity, velocity, temperature, etc.), boundary conditions,
Page 10 of 21
relaxation time and lattice spacing are defined as parameters to initialize the procedure of LB model. Initial equilibrium distribution functions are calculated through given density, velocity and temperature according to Eqs. (3) and (14). Initial distribution functions are given by f i x , t 0 = f i ( e q ) x , t 0 and Ti x , t 0 = Ti ( eq ) x , t 0 . Then the distribution functions are updated by the evolutions of velocity field and temperature field. In the freezing stage, the initial temperature of droplet is set to 0 °C, because the nucleation will release latent heat in the recalescence stage, which makes the temperature of droplet increase up to the equilibrium freezing temperature 0 °C (Hindmarsh et al., 2003) at the end of the recalescence stage. The initial volume fraction value of solid phase (fs) is set to zero based on the assumption that the initial droplet is totally liquid. At this point, the phase transition model is added and the phase of node changes from liquid to solid successively with temperature dropping. By Eqs. (1) and (13) we obtain the distribution functions in velocity and temperature field, and then the density , velocity and temperature in node can be calculated by Eqs. (5), (6) and (15). For a detailed flowchart, see Fig. 2.
2.6. Validation of LB model First the LB model will be verified elementarily by comparison with other similar studies. Since the solidification of suspended droplet has been studied and the phase transition process is similar to that of water droplet on cold flat plate, we apply the LB model to simulate the temperature variation of suspended spherical droplet and compare the data from reference. In this case, the temperature of ambient air Ta = −20.0 °C, the velocity of air v = 0.42 m/s and the droplet radius R = 0.63 mm. The simulation by LB model starts when the droplet temperature is -5.0 °C after freezing, until it cools to within 1.0 °C of the ambient air temperature. As shown in Fig. 3, the temperature curve predicted by LB model agrees well with both the published simulation results and experimental data (Hindmarsh et al., 2003). The relative mean error is 1.11 % between the present simulation result and experimental data. It is verified that the LB model and the procedure used in simulation are correct. Next, we
Page 11 of 21
will simulate the solidification process of water droplet on cold flat plate by the present LB model, and the validity of LB model will be further proved through the comparisons with our experimental results.
3. Experimental methods
To verify the present lattice Boltzmann model by experimental results, an experimental system is set up which consists of a list of components as shown in Fig. 4. The test section is placed into an isolated organic glass cavity to avoid the interference of environment, which configures two quartz glass viewing ports on both top and side of the cavity to observe droplet by the microscopic image system. The cold plate is made of red copper sized at 60 mm × 40 mm × 8 mm. The sessile deionized water droplet is distributed by a microliter syringe on the center of plate. The microliter syringe has a droplet size deviation of 0.01 μl. The semiconductor chilling plate (C-1208) is used as the chiller for the copper plate. The circulated cooling water piped from the low-temperature thermostat bath (DC-4015, ±0.1 °C) is used to cool the hot side of the semiconductor chilling plate. Contact digital temperature and humidity sensor (JWSH-5VBDD, ±0.5 °C & ± 3% RH) is installed for monitoring and recording data on humid air condition. The cold plate temperature is measured by T-type thermocouples (GBTS200, ±0.1 °C) located at five different points inside the copper plate. The exact locations of thermocouples inserted in the cold flat plate are shown in Fig. 5. All the thermocouples are located at the same level of 1 mm beneath the surface of cold flat plate. Relative mean deviations (RMD) are listed in Table 1, which indicates that relative mean deviations of the temperatures among these points are less than 1.0 %. Therefore the plate temperatures used in the boundary of LB model are obtained by averaging the measuring values of all T-type thermocouples. On the other hand, another T-type thermocouple is inserted into the droplet which can capture center temperature of
Page 12 of 21
droplet (Hindmarsh et al., 2003). The enhancement of heat transfer by thermocouple will lead to the higher measured temperature of droplet than its actual value. To correct the measured temperature, we insert thermocouple into ice-water mixture in a refrigerator and the insertion depth is equal to that in droplet. It is found that the mean temperature measured by thermocouple is about 0.45 °C higher than ice-water (0 °C). Thus the present measured temperatures of droplet are corrected according the above error analysis. The accuracy of center temperature measurement is about ±0.1 °C. All data collected in present experimental system are recorded by a data acquisition system (Agilent-34970A) to facilitate the following research.
Table 1 - Experimental data by five T-type thermocouples Time (s)
Point 1 (°C)
Point 2 (°C)
Point 3 (°C)
Point 4 (°C)
Point 5 (°C)
RMD (%)
0.0 1.0 2.0 3.0 4.0
-14.955 -15.102 -15.225 -15.336 -15.406
-14.916 -15.036 -15.149 -15.250 -15.325
-15.083 -15.216 -15.334 -15.444 -15.522
-14.806 -14.915 -14.993 -15.115 -15.196
-14.698 -14.747 -14.861 -14.998 -15.145
0.750 0.918 0.981 0.904 0.774
The microscopic image system consists of three parts: a stereomicroscope (SZX-ZB7, OLYMPUS) whose magnification is from 12 to 84, a supersensitive high-speed recording camera (MotionBLITZ EoSens mini2) which can automatically capture momentary freezing phenomenon of water droplet, and a luminescence house/annular illuminator Unit (MLC-150C, MOTIC) which can emit luminescence light without thermal radiation for observation.
4. Results and discussion
4.1. Temperature distribution variation
Fig. 6 shows the comparisons of simulated center temperatures of droplet with experimental results change over time at different cold plate temperatures. The
Page 13 of 21
temperature of air Ta = 15.0 °C, the relative humidity of air RH = 40.0 %, the contact angle θ = 84.6 °, and the volume of water droplet V = 15 μl. For three cases of plate temperatures Tw1 = −29.5 °C, Tw2 = −15.6 °C and Tw3 = −9.5 °C, the relative mean errors are 6.05 %, 4.73 % and 3.61 %, respectively. It can be seen that comparisons show acceptable agreement. According to the Fig. 6, it shows that there are four stages as described as reference (Hindmarsh et al., 2003): pre-cooling or supercooling, recalescence, freezing, and solid cooling or tempering. In the first stage, the droplet temperature keeps declining to a value below the equilibrium freezing point, which reduces faster with lower plate temperature. When the temperature is low enough for spontaneous crystal nucleation, an abrupt temperature occurs in the second stage due to latent heating. Since the growth of solid phase occurs instantaneously, the third stage is too short to be observed. Finally, the temperature continues to drop until close to that of plate. From Fig. 6, it is clearly observed that the temperature of cold plate gets lower, the freezing occurs earlier. This is because the colder plate leads to the faster thermal transport from droplet.
Fig. 7 shows the simulated temperature variations on the center axis of droplet at various time points. In the pre-cooling or supercooling stage (t < 0), the water droplet is cooled by cold flat plate. The temperature curves show a liner increase with height due to heat conduction in liquid droplet. In the freezing stage (t > 0), we can see that temperature curves show different changing trends due to phase transition. When at 1.0 s and 4.0 s, the freezing is still incomplete. A liner temperature increase along with height is observed in solid phase area. But in liquid phase area, temperature changes more gently and is influenced by ambient air close to top of droplet. When at 10.0 s, the freezing is basically completed, and therefore the temperature shows a liner increase again with height in frozen droplet.
Page 14 of 21
Additionally, the temperature distributions at various times points are plotted in Fig. 8 by Tecplot 360. We can observe the temperature field evolution of water droplet and ambient air. In Figs. 8 (a) and (b), water droplet is in the pre-cooling or supercooling stage, and the cold flat plate gradually reduces the temperature of liquid droplet to below 0 °C. In the recalescence stage, the temperature of droplet increases up to the equilibrium freezing temperature 0 °C with latent heat releasing. In Fig. 8 (c), the droplet starts to freeze (t = 0). The air temperature near the plate has already fallen to below zero, and cold air will take heat away from surface of droplet. Therefore the temperature distributions in droplet are shown as many concave isotherms and the solid-liquid interface as the one isotherm surely becomes concave with small curvature (see Fig. 10).
4.2. Volume fraction variation of solid phase
Fig. 9 shows the comparisons of simulated volume fractions of solid phase change over time at different plate temperatures with experimental data in the freezing period. In the experiment, for convenient analyzing the variation of solid phase, we assume that the water droplet maintains the shape of a spherical segment and the volume fractions of ice are obtained via measurement of height and width of ice in droplet. In Fig. 9, all three predict curves that are slightly larger than experimental data. This is because the solid-liquid interface in droplet is a concave surface with small curvature (as presented in Fig. 8) rather than horizontal surface and thus the experimental data is underestimated from measured height and width of ice. For three cases, the relative mean errors are 2.87 %, 3.00 % and 2.50 %, respectively. With the water droplet freezing, the volume fraction of ice increases approximate linearly over time at the earlier period, but the increasing rate at the end of phase transition period is slower than that in the beginning. Note that the quicker droplet freezing occurs on the colder plate. When at about 5 seconds the volume fractions of ice all achieve more than 50 % on three cases, and the freezing processes in the experiment are finished at approximately 12.62, 14.75 and 23.20 s, respectively. In Fig. 9, faster phase transition at the earlier period is due to the lower thermal resistance of cold flat plate, but over time the growing thickness of ice will lead to the higher thermal resistance which slows down the freezing rate.
Page 15 of 21
Fig. 10 shows the visualization comparisons of water droplet evolution on cold flat plate between simulations and photographs. Although from Fig. 9 we can deduce that the general shape of ice can be predicted by LB model in form of the volume fractions, more detailed comparisons of ice shape can be found in Fig. 10. As shown in simulation results, the yellow and blue area in color contour represent solid and liquid phase, respectively, and the multicolor area symbolizes status of solid-liquid mixture. In Fig. 10, compared to the freezing droplet in photographs the concave solid-liquid interfaces are also well displayed in simulation results. Furthermore, the status of solid-liquid mixture is more clearly observed in simulation results. Note that the volume of the droplet varies due to solidification in photographs. This is because more ice forms with freezing and ice has a lower density than water. Thus, the mean density of total droplet decreases with the occupancy of ice increasing, which lead to volume increase and slight deformation of droplet. In summary, via the comparisons and analyses in Figs. 6, 9 and 10, the present LB model is verified to effectively simulate the phase transition process of water droplet.
5. Preliminary study on frost formation process by LBM The LB model has been verified to effectively simulate the solidification process of water droplet. We apply the LB model to further study the frost formation process. In physical model of frost formation, the crystal nucleus forms and frost crystal grows when humid air passes above a cold flat plate having a temperature below the freezing point of water. The frost layer composed of a large number of frost crystals behaves as a porous medium into which diffusion of water vapor leads to an increase of both thickness and density (Hermes et al., 2009). In present LB model, the frost layer is considered as porous media and described with 800 × 600 lattice units. At the start, the original nuclei randomly distribute at the bottom of the domain, and then the frost crystals grow from these nuclei to the frost layer. In each time step, we update frost layer mass by nucleation theory and obtain the density distribution of frost layer. Fig. 11 shows the change of frost heights with time at different cold plate temperatures. The temperature of air Ta = 10.8 °C, the relative humidity of air RH =
Page 16 of 21
41.2 %, and the velocity of air v = 0.16 m/s. The simulation results are compared with our experimental results, and the data curves match well. It can be seen that the colder plate has the higher frost layer at the same time. When the cold plate temperature is lower, the super saturation degree and nucleation rate are higher, and thus more frost crystals are generated. But with the increase of frost height, the nucleation rate gets lower due to the rising temperature, so the growth rate of frost becomes slower gradually.
Fig. 12 illustrates the density distribution of frost layer at different time points. In this case, the cold plate temperature Tw2 = −16.5 °C. In the early-stage, the frost dendrites are formed on the cold flat plate (see Fig. 12(a)). As the frost layer grows, the bottom of frost layer becomes denser (see Fig. 12 (b) and (c)). Finally, the deeper frost crystals turn into a solid mass of ice (see Fig. 12 (d)). The densities along the height of frost layer from four time points are shown in Fig. 13. We can see that the upper frost layer has the same trend of density variation at four time points. When at 300 min frost layer growth period complete basically, the density curve becomes gentle at the bottom of frost layer and the porous structures are mainly distributed at depths of less than 1.07 mm (density < 800 kg/m3) with the porosity from 13 % to 99.7 %. After that, the frost layer grows very slowly due to increased temperature on the top. Meanwhile humid air diffuses and deposits in lower frost crystals and leads to a slow-growing ice layer.
6. Conclusion This paper presents a two-dimensional lattice Boltzmann model which is able to compute the solid phase evolution and thermal transport in the solidification of water droplet on cold flat plate. Results indicate that (1) the center temperature of droplet drops with decreasing plate temperature, and lower plate temperature (-29.5 °C compared with -15.6 °C and -9.5°C) leads to earlier occurrence of freezing, (2) the center axis temperature of droplet in the solid phase area increases linearly along with height, but temperature changes more gently in liquid phase area, (3) the temperature
Page 17 of 21
field evolution of water droplet and ambient air can be observed by the temperature cloud pictures from simulation results, (4) faster phase transition occurs at the earlier period due to the lower thermal resistance and the growing thickness of ice leads to the higher thermal resistance which slows down the freezing rate, when at about 5 seconds the volume fractions of ice all achieve more than 50 % on three cases, and (5) the shape of solid phase in freezing water droplet can be well predicted by present LB model. In addition, the preliminary study on the frost formation process indicates that (1) the colder plate (-28.6 °C compared with -22.5 °C, -16.6 °C and 8.4°C) has the higher frost layer at the same time, and (2) the growth rate of frost gets lower with the rising of frost height, frost layer growth period is complete basically at 300 min when the cold plate temperature is −16.5 °C, the density curve becomes gentle at the bottom of frost layer and the porous structures are mainly distributed at depths of less than 1.07 mm (density < 800 kg/m3) with the porosity from 13 % to 99.7 %. By comparison with our experimental data, the relative mean errors for three cases are 6.05 %, 4.73 % and 3.61 % on center temperatures of droplet and 2.87 %, 3.00 % and 2.50 % on volume fractions of solid phase, thus good agreement is observed. It is expected that the LB model is able to simulate the physical characteristics in the solidification of water droplet and establish a relationship with frosting research.
Acknowledgments The present paper is financially supported by the National Natural Science Foundation of China (No. 51106013, 50806059), and China Postdoctoral Science Foundation (No. 2012M511998).
REFERENCE Aidun, C.K., Clausen, J.R., 2010. Lattice-Boltzmann method for complex flows. Annu. Rev. Fluid Mech. 42, 439-472. Bahadur, V., Mishchenko, L., Hatton, B., Taylor, J.A., Aizenberg, J., Krupenkin, T., 2011. Predictive model for ice formation on superhydrophobic surfaces. Langmuir 27, 14143-14150. Cantor, B., O'Reilly, K., 2002. Solidification and casting. Institute of Physics
Page 18 of 21
Publishing, Bristol. Chen, S.Y., Doolen, G.D., 1998. Lattice Boltzmann method for fluid flows. Annu. Rev. Fluid Mech. 30, 329-364. Esen, M., 2000. Thermal performance of a solar-aided latent heat store used for space heating by heat pump. Sol. Energy 69, 15-25. Esen, M., Ayhan, T., 1996. Development of a model compatible with solar assisted cylindrical energy storage tank and variation of stored energy with time for different phase change materials. Energy Convers. Manage. 37, 1775-1785. Esen, M., Durmuş, A., Durmuş, A., 1998. Geometric design of solar-aided latent heat store depending on various parameters and phase change materials. Sol. Energy 62, 19-28. Eshraghi, M., Felicelli, S.D., 2012. An implicit lattice Boltzmann model for heat conduction with phase change. Int. J. Heat Mass Transfer 55, 2420-2428. Guo, Z., Shi, B., Zheng, C., 2002. A coupled lattice BGK model for Boussinesq equations. Int. J. Numer. Methods Fluids 39, 325-342. Guo, Z., Zhao, T.S., 2005. A lattice Boltzmann model for convection heat transfer in porous media. Numer. Heat Transfer, Part B-Fund. 47, 157-177. Guo. Z., Zheng, C., 2009. Theory and applications of lattice Boltzmann method, Science Press, Beijing. Hayashi, Y., Aoki, A., Adachi, S., Hori, K., 1977. Study of frost properties correlating with frost formation types. J. Heat Transfer 99, 239-245. He, Y., Wang, Y., Li, Q., 2009. Lattice Boltzmann method: theory and applications. Science Press, Beijing. Hermes, C.J.L., Piucco, R.O., Barbosa, J.R., Melo, C., 2009. A study of frost growth and densification on flat surfaces. Exp. Therm. Fluid Sci. 33, 371-379. Hindmarsh, J.P., Russell, A.B., Chen, X.D., 2003. Experimental and numerical analysis of the temperature transition of a suspended freezing water droplet. Int. J. Heat Mass Transfer 46, 1199-1213. Hoke, J.L., Georgiadis, J.G., Jacobi, A.M., 2000. The interaction between the substrate and frost layer through condensate distribution. Air Conditioning and Refrigeration Center. College of Engineering. University of Illinois at Urbana-Champaign. Jiaung, W.S., Ho, J.R., Kuo, C.P., 2001. Lattice Boltzmann method for the heat conduction problem with phase change. Numer. Heat Transfer, Part B-Fund. 39,
Page 19 of 21
167-187. Lee, K.S., Kim, W.S., Lee, T.H., 1997. A one-dimensional model for frost formation on a cold flat surface. Int. J. Heat Mass Transfer 40, 4359-4365. Martys, N.S., Chen, H., 1996. Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method. Phys. Rev: E. 53, 743-749. Matsumoto, K., Honda, M., Ito, Y., Shirai, D., 2014. Study on measurement of frost dimensions/distribution and frost crystals scraping force using scanning probe microscope (investigation on influence of humidity). Int. J. Refrigeration 38, 341-351. Miller, W., 2001. The lattice Boltzmann method: a new tool for numerical simulation of the interaction of growth kinetics and melt flow. J. Cryst. Growth 230, 263-269. Mortensen, A., Jin, I., 1992. Solidification processing of metal matrix composites. Int. Mater. Rev. 37, 101-128. Semma, E., Ganaoui, M.El., Bennacer, R., Mohamad, A.A., 2008. Investigation of flows in solidification by using the lattice Boltzmann method. Int. J. Therm. Sci. 47, 201-208. Shan, X., Chen H., 1993. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev: E. 47, 1815-1820. Shan, X., Chen H., 1994. Simulation of non-ideal gases and liquid-gas phase transitions by the lattice Boltzmann equation. Phys. Rev: E. 49, 2941-2948. Sheikholeslami, M., Gorji-Bandpy, M., Ganji, D.D., 2014. Lattice Boltzmann method for MHD natural convection heat transfer using nanofluid. Powder Technol. 254, 82-93. Tabakova, S., Feuillebois, F., 2004. On the solidification of a supercooled liquid droplet lying on a surface. J. Colloid Interface Sci. 272, 225-234. Wu, X., Dai, W., Xu, W., Tang, L., 2007. Mesoscale investigation of frost formation on a cold surface. Exp. Therm Fluid Sci. 31, 1043-1048. Yang, D.K., Lee, K.S., Cha, D.J., 2006. Frost formation on a cold surface under turbulent flow. Int. J. Refrigeration 29, 164-169. Yu, D., Mei, R., Luo, L.S., Shyy, W., 2003. Viscous flow computations with the method of lattice Boltzmann equation. Prog. Aerosp. Sci. 39, 329-367. Yun, R., Kim, Y., Min, M., 2002. Modeling of frost growth and frost properties with
Page 20 of 21
airflow over a flat plate. Int. J. Refrigeration 25, 362-371. Fig. 1 - Schematic illustration of D2Q9 (Two-dimensional, nine-speed lattices) Fig. 2 - Flow chart of the present LB simulation Fig. 3 - Comparison of center temperature variations over time Fig. 4 - Schematic representation of experimental system Fig. 5 - Locations of thermocouples inserted in the cold flat plate
Fig. 6 - Center temperature variations of droplet at different plate temperatures Fig. 7 - Simulated temperatures on the center axis of droplet at various time points Fig. 8 - Temperature contours of water droplet at various time points (Tw1 = −29.5 °C) Fig. 9 - Volume fractions of ice change over time at different plate temperatures Fig. 10 - Comparison of the water droplet evolution between simulations and photographs (Tw1 = −29.5 °C) Fig. 11 - Frost heights with time at different cold plate temperatures Fig. 12 - Density distribution of frost layer at different time points
Fig. 13 - Density variation along the height of frost layer at different time points
Page 21 of 21