A Numerical Investigation of Frost Growth on Cold Surfaces Based on the Lattice Boltzmann Method

A numerical investigation of frost growth on a cold flat surface was presented based on two-dimensional Lattice Boltzmann model. This model has been validated to have less prediction error by past experiments. According to the results, it is shown that average frost density appears different at an increasing rate at different frosting stages. In addition, cold surface temperature has great influence on frost growth parameters such as frost crystal deposition mass, frost deposition rate, and frost crystal volume fraction. It was found that the frost crystal deposition mass, frost crystal volume, and the deposition rate first increase rapidly, then gradually slow down, finally remaining unchanged while the cold surface temperature decreases. The further away from the cold surface, the more sparser the frost layer structure becomes due to the smaller frost crystal volume fraction.


Introduction
Heat exchangers will be subjected to deposition and progressive growth of frost, when humid air flows through the heat exchangers' surface where the temperature is less than both zero degrees centigrade and the ambient air dew point temperature.Because of the frost layer insulation and its blockage of airflow, frosting will bring about the reduction of a heat exchanger's performance and a greater pressure drop.In addition, unavoidable defrosting has to consume extra energy.Therefore, to guarantee normal operation of the heat exchanger, a frosting investigation on a cold surface is necessary.Apart from previous experimental studies [1][2][3], theoretical studies on frost growth predictions have been paying more attention.
A lot of predictions of frost growth have been carried out.In 1984, a molecules diffusion model was proposed by O'Neal et al. [4].In their model, the frost layer was regarded as a porous media and the mass transfer in the frost formation process not only increases frost thickness but also increases frost density simultaneously.Later, this model was continually simplified and improved by researchers [5][6][7].
Tao et al. [8] assumed the frost a homogeneous porous medium and proposed a model to simulate the fully developed frost growth defined by Hayashi et al. [9]. Lee et al. [10] assumed humid air at the frost surface as a saturated state and developed a one-dimensional model considering water vapor molecular diffusion.Different from previous theoretical models with the assumption of saturated humid air, a new model within the threshold (15%) was established by Webb et al. [11] to predict frost deposition and growth.According their work, the air at the frost surface was confirmed to be in supersaturated state.Hermes et al. [12] formulated a mathematical model within ±10% error bands.In their model, they considered the air to be in a supersaturated state and regarded the frost as a porous media.In another study, a one-dimensional frost growth model was developed by Kandula [13] based on Jones and Parker [14] and Cheng et al. [15].New frost correlations for density and effective thermal conductivity on a flat plate were formulated in Kandula's study.Depending on semi-empirical correlations, Loyola et al. [16] presented a one-dimensional model to simulate the frost formation on parallel-plate channels under the conditions of the supersaturated humid air near the frost surface.All these above listed studies are one-dimension models simulating frost growth, assuming the frost layer to be a homogeneous porous media.These models mainly depend on empirical correlations to calculate frost parameters such as density and effective thermal conductivity and so on.Furthermore, some two-dimensional models were established to study frost growth without using frost empirical correlations in recent years.Lüer et al. [17] developed a two-dimensional transient model to study frost formation between parallel plate channels based on the conservation laws of mass, energy, and momentum.Their numerical results show good agreement with their experimental findings.With this model, local temperature and porosity of frost layer can be provided.Armengol et al. [18] also presented a two-dimensional model for frost growth over parallel plates adopting the finite volume method.The predicted frost thickness shows agreement with experimental data within a 10% error.Without adopting any experimental correlations, Lee et al. [6] established a transient two-dimensional prediction of frost deposition on a cold flat surface with agreement of about 10%.Another two-dimensional model was also proposed by Lenic et al. with initial frost density of 30 kg/m 3 and initial frost thickness of 0.02 mm [19].In his model, the moist air was assumed to be in a super-saturated state.Some CFD models were also adopted to simulate frost formation within acceptable deviation [20][21][22].More recently, Negrelli et al. [23] investigated the frost effective thermal conductivity in details based on fractal theory with an approximately ±15% uncertain threshold compared to experimental results.Although the above models show good tendency with experimental points, the deviations in most models are about 10%.Therefore, better accurate prediction is needed.
The lattice-Boltzmann (LB) method is a novel numerical method.Compared with the traditional simulation methods, the LB model has attracted wide attention for its adaptability to complex structures and has achieved success in some areas, such as porous media and multiphase [24].Manz et al. [25] conducted experiments and numerical simulations of randomly stacked porous bodies and found that the results of the LB simulation agree well with the experimental data, demonstrating the suitability of LB for simulating flow in porous media.Tang et al. [26] proposed a simplified analytical model for micro-scale porous media, and used the LB model to simulate the flow field of micro-scale porous media.Li et al. [27] proceeded from the lattice-Boltzmann collision model to realize the numerical simulation of the fluid migration process in dual porous media.Xu et al. [28] adopted a regular hexagonal model (FHP-7Bit model).They also established a corresponding seepage numerical model according to the relevant LB theory.This model was verified to be accurate.Wang et al. [29] used the LB model to simulate the process of gas seepage in dense porous media.The correctness of the simulation results was verified.With this approach, the flow field in the porous medium in the seepage flow was obtained.All the above studies show that the LB model is feasible and suitable to simulate the flow and heat transfer phenomenon in porous media.As mentioned above, it is known that the frost layer, consisting of a great number of frost crystals, acts as a porous media [4].We have been trying to use the LB model to predict frost growth in recent years.We established a LB model to simulate the solidification process of water droplets effectively.Furthermore, we have also adopted this LB model to investigate frost growth and have acquired some preliminary findings [30].In order to obtain more detailed information on frost growth, we will continue our study deeply following our previous LB model.
This study presents a two-dimensional LB model to simulate the frost growth on a cold flat surface based on the fractal theory (DLA model).In this model, the frost layer was regarded as a heterogeneous medium and the humid air near the frost surface is considered to be in actual super-saturation state.Compared to experimental data, the simulation results show good agreement; the maximum relative error is 8.05%.According to this model, the dynamic variation rules of frost layer properties have been investigated thoroughly.
The frost layer growth model based on the LB method has a good agreement with the experimental results.The model can be used to predict the growth of the frost layer and explain the characteristics of each stage of frost layer growth.It can provide more detailed information about frosting on flat cold surfaces.Meanwhile, it can also help to find reasonable defrost strategies and to optimize heat exchanger design under frosting conditions.In the following research, this study will provide necessary information for frosting issues on complex cold, instead of flat cold, surfaces.

Numerical Model and Theory
This simulation was carried out based on the following hypothesis: (1) for the nucleation process of frost crystals, homogeneous nucleation is adopted; (2) the crystal nucleus with a critical size formed by the water vapor condensation directly deposited on the existing frost crystals; (3) the humid air is considered as an incompressible Newton fluid; (4) in the simulation process, super-saturation and the relative humidity of humid air are assumed to be fixed; (5) during the frost growth, the melting and sublimation of frost crystals in some areas were neglected.
In order to create frost porous structure, the fractal theory (DLA model) is adopted in this two-dimensional LB model.According to DLA model, a certain number of seed particles are randomly generated on the cold surface (X-axis) firstly.Then, in order to simulate the frost crystal after the phase transition nucleus of water vapor, a certain number of random particles are generated in the calculation region using a random function; the probability of random particles generated is related to the freezing probability in Equation (1).The distance between random particles and seed particles is firstly calculated, and random particles will attach to the main frost crystal at the seed position with the smallest distance.In this way, the random particles generated at the frost top layer calculation area will preferentially attach to the main frost trunk, which is the strong longitudinal growth characteristics.At the same time, in order to simulate the water vapor bypasses around the main frost crystals and condenses at the inside of the frost layer, random particles are also generated inside the frost layer and move towards downward, left, and right directions with equal probability.Then these random particles will attach to the nearest cold surface or existing frost crystals.

The Mass of Frost Layer
Based on the nucleation theory, the increase of frost crystal mass within the controlled volume per unit time is calculated by the following Equation [31]: where the critical radius of the nucleus r c is: where σ = (99.5 − 0.075T) × 10 −3 J/m 2 is the interface energy, α = p v /p s is the water vapor saturation ratio, S = α − 1 is the super-saturation degree of humid air.The nucleation rate J is: where k B = 1.38 × 10 −23 J/K is the Boltzmann constant.Based on the assumption of homogeneous nucleation, the contact angle θ is 180 • .The kinetic constant of nucleation rate J 0 is: where M = 2.989 × 10 −26 kg is the water molecular mass.The non-isothermal correction factor ζ is: In Equation ( 1), P = T 0 − T f s /(T 0 − T w ) is the freezing probability, T 0 = 273.15K.ε i is the frost crystal volume fraction.Based on the frost crystal volume fraction ε i , the parameter µ(ε i ) [32]  is written as The diffusion resistance factors C 0 , F µ and k are equal to 0.58, 3, and 10, respectively.The frost crystal volume fraction ε i is written as The equivalent density of frost layer is defined as porosity which refers to the volume fraction of voids between particles in a porous structure.

The Velocity Field and the Temperature Field in LB Model
In this LB model, the Bhatnagar-Gross-Krook (BGK) [33,34] approximation with two dimensions and nine velocities was adopted.For the velocity field, the distribution function is: where f i is the velocity distribution function.x and t are the position of the particle and time, respectively.τ is the dimensionless relaxation time which is connected with the kinematic viscosity ν = c 2 s (τ − 0.5)δt.The discrete velocity e i is: where δx and δt are the lattice spacing and the time step, respectively.The equilibrium distribution function of velocity field f (eq) i is: where c s = (δx/δt)/ √ 3 is the isothermal sound velocity.The weight coefficients ω i are {4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36}.
The macroscopic density ρ and velocity u are defined as: where v is an auxiliary velocity with two parameters c 0 = 1 + ε p δtν/2K /2 and c 1 = ε p δtF ε / 2 √ K .
F ε = 1.75/ 150ε 3 p is the geometric function and K = ε 3 p d 2 p / 150 1 − ε p 2 [35] is the permeability, d p is the diameter of the solid particle.
For the temperature field, the distribution function is: where the dimensionless relaxation time τ c is connected with the thermal diffusivity α m = c 2 s (τ c − 0.5)δt.The discretized source term is φ t = ω i h sg ε i − ε i /c ps δt.The equilibrium distribution function of temperature field T (eq) i is: The macroscopic temperature T is: where the heat capacity ratio [36] ϑ = ε p + 1 − ε p ρ s c ps /ρ v c pv reflects the influence of the porous medium on the temperature field.ρ s and ρ v are the solid and fluid densities, c ps and c pv are the solid and fluid specific heats at the constant pressure.
The evolution of LB model mainly includes two processes: migration and collision.

1.
Migration: Collision: During the process of evolution, temperature will be transferred into the system through the migration of particles on the boundary.And particle collision will occur due to the change of temperature.Therefore, migration process is regarded to happen before collision in this model.

Initial and Boundary Condition
The initial condition.The initial velocity of the velocity entrance boundary is equal to the air velocity, and the other nodes are zero.The initial temperature of cold wall is equal to the cold surface temperature, and the other nodes are equal to the air temperature.
The boundary condition.In the rectangular calculation area, the macroscopic temperature and velocity on the corresponding boundary were set as following.The velocity of left boundary is constant and equal to the inlet air velocity, and the temperature is equal to the humid air temperature.The velocity and temperature of the nodes on the right boundary and the top boundary are determined by the nearby front layer nodes by the non-equilibrium extrapolation format in the calculation process.The velocity of the bottom boundary (cold surface) is zero and the temperature is constant and equal to the cold surface temperature.
In the LB model, the velocity and temperature of the boundary nodes in each direction are determined by the distribution functions.Based on the macroscopic parameters of the boundary, the boundary conditions of the corresponding nodes are obtained.The details are as follows.For the cold wall nodes, the bounce-back format was adopted.A detailed schematic is shown in Figure 1a, and the distribution function of the boundary node (i, 1) is: For the other three boundary nodes, the non-equilibrium extrapolation format was adopted.A detailed schematic is shown in Figure 1b, the nodes E, O, and A are located on the boundary, the nodes B, C, and D are located in the flow field, the nodes F, G, and H are located outside the flow field.The distribution function of the boundary node O is calculated in the following Equations:   In the simulation process, the initial physical parameters were firstly determined.The physical parameters include humidity, velocity, temperature and density of the humid air, cold surface temperature, relaxation time, lattice space, and boundary conditions.At the same time, a certain number of random nuclei were distributed on the cold wall nodes.Then the particles representing the frost crystal after the phase transition nucleus of water vapor, are randomly generated at the top of the calculation area and these particles move towards the left, right, and bottom directions with equal probability until they contact with the solid nodes.According to the frost crystal volume fraction, the number of crystals can be calculated.The initial velocity and temperature distribution functions are determined by f i (x, t 0 ) = f i (eq) (x, t 0 ) and T i (x, t 0 ) = T i (eq) (x, t 0 ), respectively.Based on the Equation ( 1), the mass increase of the frost layer can be calculated.At the same time, the frost crystal volume fraction will be updated.Finally, based on the boundary conditions, the velocity and temperature distribution functions, as well as the macroscopic parameters, such as density, velocity and temperature, can be obtained.

Validation of LB Model
Figure 3 shows the comparison of the model results with the experimental data provided by a previous paper [19].In Figure 3, air temperature T a = 19.8• C, air relative humidity ϕ = 58.0%,air velocity u = 0.6 m/s, and cold surface temperature T w = −20.5 • C. It can be seen from Figure 3 that the simulation results are in good agreement with the experimental data and that the maximum relative error is 8.05%.

Results and Discussion
Here, a domain of 4 mm × 3 mm is used, and the lattice number adopts 800 × 600.The following simulated air conditions are set as followings: T a = 10.8 • C, ϕ = 41.2%, and u = 0.16 m/s.

Frost Density
Figure 4 shows the variation of average density of frost layer over time.The cold surface temperature T w = −16.5 • C. Generally speaking, the average density of frost layer increases with time.However, it can be seen that frost average density shows a different increasing rate at different frosting stages.In the incipient frosting, the curve has a little concave tendency.It is known that the frost exhibits strong dendritic growth during this early stage of frosting.That means the water vapor from the humid air mainly contributes to increase frost thickness while the rest proportion of water vapor increases and the density decreases.Therefore, the frost density growth slows down in the early stage of frosting.As the frost crystals continue to grow, the water vapor molecules continue to diffuse and desublimate to form nucleation deposits inside the frost layer.Therefore, the average density of frost increases gradually over time.However, due to the rise of the frost crystal surface temperature at the later stage of the frost growth, the freezing probability of the crystal nuclei decreases.As a result, the increase of the frost layer average density in the later stage is slower.It can be seen that the frost thickness increases with the decrease of cold surface temperature and that the lower the cold surface temperature is, the greater the frost thickness increases.The main reason is that the water vapor will undergo the phase transition and form the crystal nucleus when the humid air approached the cold surface.Therefore, there will be a water vapor concentration gradient near the cold surface, which can force water vapor to move toward the cold surface.The lower the cold surface temperature, the greater the driving force due to the concentration gradient.It means that more water vapor molecules are able to contact and deposit on the cold surface or the existing frost crystal.Consequently, frost grows dramatically faster with the decrease of cold surface temperature.

Frost Deposition
The quantitative variation of frost crystal deposition can be used to measure the frost growth speed.Figures 6 and 7 show the effect of cold surface temperature on the frost crystal deposition and the frost crystal deposition rate.Besides, it can also be found that the frost crystal deposition rate increases with time and finally tends to remain unchanged.During the early frosting period, the increase in the frost crystal deposition rate is very rapid and occurs in a short time.This is due to the fact that in the early frost period, the crystals continue to grow upward in a dendritic form.Such dendritic structure acts as fins enhancing the effect of heat transfer and increasing the probability of phase transition.As a result, the frost crystals deposit faster.At the later frosting period, the water vapor hardly diffuses into the interior of the frost layer due to the increase of the frost surface temperature.Such a diffusion opportunity lessens with frost developing.Therefore, it is reflected that the deposition rate of frost crystals slows down gradually and approaches a constant value.Therefore, the deposition rate of the frost crystal is slowed down.

Frost Crystal Volume Fraction
Frost crystal volume fraction is the volume percentage of the frost crystal in the frosted area.It directly reflects the number of frost crystals and the density of the frost layer.At the same time, it can also reflect the porosity of the frost structure.The closer to the cold surface, the larger the frost crystal volume fraction, the more concentrated the frost crystal, and the larger the frost structure density.The further away from the cold surface, the smaller the frost crystal volume fraction, and the more sparse the frost crystal distribution.As the frost thickness increases, the volume fraction gradually decreases.This is due to the fact that water vapor diffuses from the outside to the inside of the frost layer and finally deposits at the bottom of the frost layer.So the more upwards along the frost thickness direction, the more sparsely the distribution of frost crystal is.At the same frost thickness, the lower the cold surface temperature, the frost crystal volume fraction is larger.

Conclusions
This paper presented a two-dimensional LB model to simulate frosting process on a cold surface based on the fractal theory (DLA model).This model was validated by experimental data with greater accuracy.According to the model, the dynamic variation rules of frost properties have been studied in detail.Some conclusions have been drawn out.
(1) Frost average density shows different increasing rates at different frosting stages.The increase in the frost layer average density in the later frost growth stage is slower due to the rise of frost surface temperature.(2) Due to the effect of the cold surface temperature on the water vapor concentration gradient, the frost thickness increases dramatically with the decrease in the cold surface temperature.(3) The frost crystal deposition mass and deposition rate are greatly affected by the cold surface temperature.With the decrease in the cold surface temperature, the frost crystal deposition mass increase rapidly, and the deposition rate first increases rapidly, then gradually slows down, finally remaining unchanged.(4) The frost crystal volume fraction is greatly affected by the cold surface temperature.When the cold surface temperature decreases, the frost crystal volume fraction is greater.The further away from the cold surface, the smaller the frost crystal volume fraction.As a result, the frost layer porous structure appears more obvious.

Funding:
The present paper is financially supported by the National Natural Science Foundation of China (No. 51776149, 51106013), which is gratefully acknowledged.

Conflicts of Interest:
The authors declare no conflicts of interest.

Figure 2
Figure 2 shows the calculation flowchart.

Figure 2 .
Figure 2. The flow chart of the Lattice-Boltzmann method.

Figure 4 .
Figure 4. Variation of average density of the frost layer over time.

Figure 5
Figure 5 shows that the frost thickness changes with time under different cold surface temperatures.The cold surface temperature T w are set as: 4.0 • C, −8.0 • C, −12.0 • C, and −16.0 • C.It can be seen that the frost thickness increases with the decrease of cold surface temperature and that the lower the cold surface temperature is, the greater the frost thickness increases.The main reason is that the water vapor will undergo the phase transition and form the crystal nucleus when the humid air approached the cold surface.Therefore, there will be a water vapor concentration gradient near the cold surface, which can force water vapor to move toward the cold surface.The lower the cold surface temperature, the greater the driving force due to the concentration gradient.It means that more water vapor molecules are able to contact and deposit on the cold surface or the existing frost crystal.Consequently, frost grows dramatically faster with the decrease of cold surface temperature.

Figure 5 .
Figure 5. Variation of frost thickness over time.

6
shows the evolution of frost crystal deposition per unit area over time at different cold surface temperatures.The cold surface temperature T w are set as: −5.0 • C, −10.0 • C, −15.0 • C, and −20.0 • C. It can be observed that as the temperature of the cold surface decreases, the frost growth rate tends to increase, and the mass of the frost crystals increases significantly.

Figure 6 .
Figure 6.Variation of frost deposition over time at different cold surface temperatures.

Figure 7
Figure 7 shows the evolution of frost crystal deposition rate with time.The cold surface temperature T w are set as: −5.0 • C, −10.0 • C, −15.0 • C, and −20.0 • C. At first sight, the cold surface temperature has a prominent effect on frost crystal deposition rate.The frost crystal deposition rate increases dramatically with cold surface temperature decreasing.Besides, it can also be found that the frost crystal deposition rate increases with time and finally tends to remain unchanged.During the early frosting period, the increase in the frost crystal deposition rate is very rapid and occurs in a short time.This is due to the fact that in the early frost period, the crystals continue to grow upward in a dendritic form.Such dendritic structure acts as fins enhancing the effect of heat transfer and increasing the probability of phase transition.As a result, the frost crystals deposit faster.At the later frosting period, the water vapor hardly diffuses into the interior of the frost layer due to the increase of the frost surface temperature.Such a diffusion opportunity lessens with frost developing.Therefore, it is reflected that the deposition rate of frost crystals slows down gradually and approaches a constant value.Therefore, the deposition rate of the frost crystal is slowed down.

Figure 7 .
Figure 7. Variation of frost crystal deposition rate over time at different cold surface temperatures.

Figure 8
shows the distribution of the frost crystal volume fraction versus different frost thickness under different cold surface temperatures.The cold surface temperature T w are set as: −5.0 • C, −10.0 • C, −15.0 • C, and −20.0 • C. It can be observed the frost crystal volume fraction distributes non-uniformly within the frost layer.

Figure 8 .
Figure 8. Distribution of frost crystal volume fraction in the frost layer.