，，2，，2，，2*
1.College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，P.R.China；
2.State Key Laboratory of Mechanics and Control of Mechanical Structures，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，P.R.China
Abstract: The impact and freezing of micro-sized droplets on cold surface is simulated by the developed numerical methods which couple the multiphase lattice Boltzmann flux solver to simulate the flow field，the phase field method to track the droplet-air interface，and the enthalpy model to determine the liquid-ice interface.The accuracy and reliability of the numerical method are validated by the comparison between the predicted morphology of the droplet impact and freezing on the surface and that from the experiment.The dynamic freezing process is investigated considering the effects of the droplet size，the impact velocity and the temperature of the cold surface.The results show that the freezing of the droplet bottom inhibits the rebound after the droplet spreading，and it may even form a hat-like shape.For the droplet with higher velocity，the ice develops faster in the radial direction and the heat transfer between the droplet and surface is enhanced.In addition，the temperature governs the dynamic behavior of the droplet center.When the surface is colder，it may form a crater in the center.The analysis on the temperature distribution inside the droplet shows that the heat flux decreases with the increasing distance to the cold surface.Moreover，with the ice growing，the decreased temperature in symmetric axis is not proportional to the surface temperature.The dimensionless temperature inside the ice becomes lower for the colder surface due to the increased temperature difference.
Key words：aircraft icing；
supercooled droplet；
impact；
freezing；
numerical simulation
When the aircraft passes through the region with micro-sized supercooled droplets during takeoff and landing，the ice accretion may occur on the airframe［1-2］.The aircraft wings，tail，engine and other instruments for detecting the meteorological parameters of the flight will suffer from the icing threat，which may deteriorate the aerodynamic performance，i.e.the decrease of the thrust and the increase of the drag［3-4］.Besides small droplets，largersized droplets are paid much attention by researchers.The dynamic behaviors of micro-sized droplets and freezing on cold surface are complicated due to the heat and mass transfer as well as the dynamics of droplets，ice and air.The study on the momentum and energy of the droplet impact and freezing helps to understand the underlying mechanism［5］.Furthermore，from microscopic level to macroscopic level，the investigations on the droplet freezing could provide reference for devising efficient instruments of anti-icing and defrosting［6］.
The impact of droplets and freezing on cold surface involves two aspects，the dynamics of droplets and heat transfer between the droplet and surface and that between the liquid and ice.The previous studies mainly focused on the individual process，and the interactions of the momentum and heat transfer were investigated by several researchers［7-10］.
When the droplet impacts on the surface，it may experience spreading，retraction，rebound or oscillation stage［11-12］.Some experimental work and simulations were implemented and the dynamic behaviors were analyzed.Li et al.［11］simulated the droplet impact onto the surface with the level set method and suggested that the droplet spreads faster on the surface when the impact velocity is higher，and the maximum spreading factor increases with the temperature.In addition，the droplet impact splash model was proposed.Zhang et al.［13］used the direct numerical simulation to elucidate the impact forces of water droplet falling on the surface.They concluded that there are two prominent peaks，the first peak occurs due to the inertial shock while the second occurs before the droplet rebound.And according to the amplitude of the second peak，different regimes were identified.Besides this，Zhou et al.［14］used the lattice Boltzmann method to simulate the group of droplet collision on the surface and illustrated that the viscosity of the droplet and air，the surface tension and other material properties may dominate the interactions of the droplets during the impact and coalescence.In experimental work，Palacios et al.［15］studied the splashing/deposition threshold for droplet impact onto the dry surface.They reported that the viscosity promotes the splashing and thin film lifting，although it inhibits the formation of the secondary droplet from the lifted thin film rim.What’s more，Ebrahim et al.［16］focused on the impact of the droplet propelled by an air stream.It was found that the gas stream did not play an important role in the spreading.However，the droplet receding was decreased by the gas.Moreover，Guo et al.［17］investigated the oblique impact dynamics of droplets.They showed that the inclination angle of the surface has little effect on the jet emission and bubble entrapment in the droplet.And they classified the jetting phenomenon into different types，which are related to the air cavity collapse.
For the heat and mass transfer，the cooling from the surface is the fundamental reason.When the droplet begins to freeze，the droplet first experiences the nucleation stage and rapidly completes the recalescence［18-20］.The latent heat is released during the phase transition.In numerical studies，Yao et al.［21-23］simulated the freezing process by the volume of fluid and enthalpy porosity model.They indicated that when the range of the local initial ice fraction was small，the results predicted by the average and local initial ice fraction showed the similar propagations of the moving liquid-ice interface.Wang et al.［24］introduced the VOSET method which combines the volume of fluid and level set method to track the interface between the droplet and air.They stressed the importance of consideration on the supercooling effect and found that the freezing front moves slowly and the freezing speed is decreased due to the greater thermal resistance of the ice.Vu et al.［25-27］adopted the front tracking method to capture the interface and demonstrated that the freezing is strongly affected by the Stefan number，the thermal conductivity and heat capacity ratios of the ice to liquid.Increase of these parameters hastens the growth rate of the liquid-ice interface.In experimental researches，Chaudhary et al.［20］presented the experimental study on the static droplet freezing on the cold surface with different wettabilities.And they demonstrated that the time taken for droplet freezing depends on the temperature at the pre-recalescence instant.And the temperature inside the droplet increases to the freezing point when the recalescence is complete.Zhang et al.［28-29］studied the single sessile droplet freezing on the surface with focus on the height of the ice inside the droplet.They described that the droplet size governs the increasing rate of the ice height.What’s more，they proposed an empirical model to estimate the required time for the droplet to complete freezing，considering the droplet size，the surface temperature and contact angle.Besides this，Stiti et al.［30］experimentally measured the ice front within the droplet and observed the spherical shape.Two materials were tested in their experiments and it was explained that the radical difference of the freezing rate was due to the thermal diffusivities.
Despite the studies introduced above investigate the droplet impact and freezing on cold surface，the complicated momentum，heat and mass transfer are supposed to be further illustrated to provide more details on the dynamic mechanism and freezing characteristic during this impact-freezing coupled process.The aim of this paper is to develop the numerical method to better track the interface of droplet-air and determine the liquid-ice interface.In addition，to reveal the underlying mechanism，different parameters are considered to elucidate their effects on the dynamic interactions and heat transfer among the droplet，ice and air.The detailed analysis on the ice formation and ice growth helps to establish the mathematical and physical model for the aircraft icing prediction in macroscopic level and to devise new devices for enhancing or preventing icing in related applications and technologies.
The rest of this paper is organized as follows.In Section 1，the numerical methods are introduced.And the computational setup and the comparison between the computational result and experimental data are presented in Section 2.Section 3 describes the effects of droplet size，impact velocity and temperature of cold surface on the droplet dynamics and heat transfer，respectively.Finally，some main conclusions are drawn in Section 4.
The micro-sized supercooled droplet impact and freezing on cold surface includes complicated dynamic and thermal physical mechanisms.There are three phases during the whole process，i.e.，liquid water and ice as well as air.To distinguish different phases and their flow characteristics，the computation of the flow field is supposed to be coupled with determinations of the liquid-air interface and liquidice interface.The developed numerical methods are introduced below.
1.1 Governing equations of the flow field
The volume of water is expanded during the impact and freezing process of micro-sized droplet on cold surface.When the phase transition occurs，the formed ice accumulates and the momentum disappears gradually.The multiphase region occupied by the liquid and air is assumed to be incompressible and the governing equations of the flow field could be written as
wheretis the time，pthe pressure，ρthe density，uthe velocity，csthe speed of sound，μthe dynamic viscosity，andFsthe surface tension.The subscripts“l” and “s” represent the variables of liquid and ice，respectively.Cis the non-air fraction，representing the volume of the fluid，which takes 1 in non-air region and takes 0 in air region.Andγis the liquid fraction which takes 0 in the ice region and takes 1 in non-ice region.The right-hand side of Eq.（1）indicates the volume expansion of the droplet.The liquid fraction is coupled with the enthalpy equation of the multiphase flow.
1.2 Phase field method
Besides the flow field，the interface between the liquid and air is captured by the phase field equation，which is shown as
It should be noted that the right-hand side of Eq.（3）accounts for the droplet volume expansion effect which plays a role in the interface between the liquid and air.
In the phase field method，the surface tensionFsis related toC，and the following form of the surface tension is applied［31-32］
whereμCis the chemical potential and can be calculated as［33］
whereβandkfare related to the surface tensionσand the width of the interfaceWin the simulation，yieldingβ=12σ/Wandkf=3σW/2.For the wettability of the cold surface，the static contact angle model is applied，which is described as［31，34］
wherenis unit vector of the surface，Cwthe non-air fraction at the cold surface andφeqthe equilibrium contact angle.
1.3 Enthalpy equation for liquid freezing
The droplet impact on the surface could be simulated by the numerical methods introduced above.The freezing process of the micro-sized droplet requires the thermal computation.The energy equation in enthalpy form is expressed as［29，35-36］
whereHis the enthalpy of the fluid，Tthe temperature andkthe thermal conductivity.During the freezing process of the droplet，the latent heatLis released and it is incorporated in the enthalpyH，which is defined by［35］
whereCpis the specific heat capacity.Due to the difference of the multiphase material properties，ρ，Cpandkshould transit smoothly from liquid to air，ice to liquid and ice to air.Therefore，the material properties can be taken by
whereψrepresentsρ，Cpork.The subscripts “l”，“s” and “a” denote the properties of liquid，ice and air，respectively.
In addition，the physical relationship between the liquid fraction and enthalpy should be established to close Eqs.（1—9）.In fact，the supercooled droplet impacts on the cold surface and the nucleation as well as the recalescence occurs with the increase of the droplet temperature to the freezing point 0 ℃.This stage is rather rapid and the cost is the sacrifice of the latent heat of the liquid.After this stage，the droplet starts to freeze，and the temperature of the liquid begins to drop below the freezing point，which is simulated in this study.The liquid fractionγis regarded as the function of the temperatureTand expressed as follows［10，35］
wherekTis the parameter to control the temperature range of phase transition and increase the numerical stability of the computation.Here，Tshiftis used to ensure that the liquid temperature is at freezing point and matches the corresponding initial liquid fractionγ0=1-CpΔT/L，where ΔTis the supercooling of the droplet.
To ensure the frozen ice stays static state，the velocity of the ice region should be maintained at zero.With the definition of non-air fractionCand liquid fractionγ，the different phases could be identified.1-C，CγandC（1-γ）denote air，liquid and ice，respectively［37］.In the ice region，the velocity is corrected by
Obviously，whenC（1-γ）equals 1，the velocity is zero and the heat transfer in the ice is only thermal diffusion.
1.4 Numerical scheme
The simulation of the droplet impact and freezing is implemented in the cylindrical coordinate system.The governing equations of the flow field is solved by the multiphase lattice Boltzmann flux solver（MLBFS）［31-32，38］，where the finite volume method is applied and the flux through the cell interface is computed by the lattice Boltzmann method.To track of the interface of liquid and air，the phase field equation is discretized with the fifth order WENO scheme to increase the accuracy of the solution［39-40］.The energy equation is solved by the finite difference method，where the second upwind scheme is adopted.The flow variables and temperature field are updated with the third order total variation diminishing（TVD）Runge-Kutta scheme［40］.
In this section，the computational setup is introduced and the comparison between the computational result and experimental data is reported to validate the accuracy and reliability of the numerical methods.
2.1 Computation setup
In the simulation，due to the axisymmetric property of the droplet impact and freezing，only half domain（400×400 cells）is adopted to balance the accuracy and computational load.The computational setup is presented in Fig.1，where the supercooled droplet impacts the cold surface with velocity ofUand the cold surface is kept with constant temperatureTwwhile the air temperature is at freezing point which is equal to that of the droplet after the nucleation and recalescence stage.The left side of the domain is axisymmetric while the right and top sides of the domain are defined as the Neumann boundary.For the cold surface，the non-slip boundary is implemented.The droplet diameterDis represented by 200 cells.
Fig.1 Schematic of the computational domain
Additionally，the properties of the liquid，air and ice are listed in Table 1［41］.Unless otherwise specified，these parameters are adopted in the following simulations.
Table 1 Properties of liquid,air and ice[41]
2.2 Validation of the freezing model
The droplet impact and freezing on the cold surface is simulated by the developed numerical methods introduced above and the result is compared with experimental data［42］as shown in Fig.2，where temperature of the droplet is -5 ℃ while the initial droplet diameter is 2.84 mm，the impact velocity is 0.7 m/s and the cold surface is set to be -30 ℃.
Fig.2 Comparison of simulation results and experimental data[42]
This condition is suitable for validation of the present numerical method to prove the tight coupling of the mass，momentum and energy for droplet impact and freezing process，as the temperature of the cold surface and the droplet supercooling are within the range of aircraft icing condition［1］.Furthermore，the Reynolds numberRe=ρlUD/μland the Weber numberWe=ρlU2D/σfor this condition areRe=1 110.6，We=18.4，respectively，which are within the range of droplet impact condition in aircraft icing［1］.And in this condition，the droplet spreading，retraction and recoil stages are shown，which are accompanied with dynamic freezing process.In Fig.2，it is evident that the droplet first spreads on the surface and the contact area is increased.When the droplet spreads to be like a pancake，the droplet begins to retract and the top of the droplet jumps.However，due to the freezing of the droplet bottom，an ice layer inhibits the detachment of the droplet from the cold surface，which was observed in the experiment［42］，and this detail is also captured by the simulation.When the kinetic energy of the droplet dissipates gradually，the droplet becomes a pancake again.
A parameter related to the deformation of the droplet is introduced here，the spreading factor，which is defined as the ratio of the droplet spreading diameterDto the initial droplet diameterD0.The dimensionless spreading factors predicted by the simulation and experiment are described in Fig.3.Despite the small discrepancy in the maximal spreading diameter，the relative error is no more than 5%.And finally，the droplet spreading factor maintains at a constant which is in accordance with the experiment and the moving of the ice caused by the numerical error［43］disappears.The qualitative and quantitative comparison results show that the simulation result is in good agreement with the experimental data，validating the accuracy of the numerical method.
Fig.3 Spreading factor predicted by simulation and experiment from Ref.[42]
The droplet impact and freezing on cold surface is governed by several conditions，including the diameter of the micro-sized droplet，the impact velocity of the droplet，and the cooling of the cold surface.To explore the effect of these factors on the dynamic behaviour of the droplet and its heat transfer with air and ice，a series of simulations are performed and the simulation conditions are presented in Table 2.It should be noted that the supercooling of the droplet ΔTis equal to the absolute value ofTw，i.e.，ΔT=|Tw|.When the supercooled droplet impacts on the cold surface，the recalescence occurs and completes instantly while the temperature of the liquid droplet immediately rises to the freezing point 0 ℃ with the cost of the partially release of the latent heat before the freezing of the droplet.
Table 2 Simulation conditions
To investigate the general dynamic and thermal characteristics of the droplet impact and freezing on the cold surface，two important dimensionless variables are introduced here，i.e.，dimensionless timet* and dimensionless temperatureθ，which are defined as
And the following discussions on the simulation results are based on these dimensionless variables.
3.1 Effect of droplet size
In the aircraft icing conditions，the size of the droplet usually distributes from 101μm to 102μm.Different-sized droplets may present various dynamic behaviors and freezing features.The evolutions of the morphologies of different-sized droplets and ice growing processes inside droplets are shown in Fig.4.
During the spreading process of the droplet，an ice layer is formed at the bottom of the droplet.With the increase of contact area of the droplet and the cold surface，the ice layer develops not only in the vertical direction but also in the radial direction.For small-sized droplets，after the spreading stage，the droplet tends to recoil and jump out of the surface att*=2.5，as shown in Fig.4（a）.As the bottom is attached to the cold surface，the adhesion force restricts the tendency to jump.The inertial force competes with the adhesion force and gradually weakens.After long-time oscillations due to the transition of kinetic energy and surface energy，the droplet tends to be flat att*=4.0，as shown in Fig.4（a）.While for the droplet with larger diameter，the droplet spreads more because of the inertia，and the increased contact area enhances the heat transfer between the droplet bottom and the cold surface.Hence，the ice layer forms in the radial direction and inhibits the retraction of the droplet，and the liquid in the droplet center becomes more flat and even concave att*=4.0 as shown in Fig.4（c）.
Fig.4 Evolutions of morphologies of different-sized droplets and ice growing processes for Cases 1 to 4
It is noted that when the droplet diameter is 200 μm（Case 4），the spreading of the droplet is severely limited by the icing process of the droplet and the whole icing process is accelerated compared with those of smaller-sized droplets（Cases 1 to 3）.The strong heat transfer at the surface causes the rapid icing and the droplet upper part tends to spread.Actually，in the spreading stage，the inertial force promotes the momentum in the normal direction transfer to the radial direction，during which，besides the viscous dissipation，the freezing process increases the dissipation too.In dimensionless time，the freezing process competes with the spreading tendency due to the inertial force and the icing is expedited for the large droplet（Case 4）so that the bottom of the droplet is frozen rapidly first and with the increased ice thickness，the freezing is slowed down because of the decreased thermal flux density in the normal direction.Att*=0.5，the droplet starts to present the obvious characteristics of spreading like Cases 1 to 3.And finally，the droplet shape becomes like a hat with thin rims att*=1.5，as shown in Fig.4（d）.
The icing process is the result of interactions between the latent heat release and cooling from the cold surface.The interface between the ice and liquid is concave or flat，as shown in Fig.4，depending on the heat transfer at the interface.The temperature evolutions inside the droplet reflect the heat transfer between the droplet bottom and the cold surface and that between the droplet and air.In Fig.5，the dimensionless temperatureθin the symmetric axis is presented.It can be seen that the temperature inside the droplet gradually decreases with time due to the cooling from the surface.And the temperature distribution in the vertical direction is not linear.From the ice bottom to the top，the temperature first increases fast and then the acceleration declines，indicating that the heat flux decreases with the increasing distance to the cold surface.This is attributed to the increasing thickness of the ice，where the heat transfer is abated near the liquid-ice interface.
Fig.5 Distribution of the dimensionless temperature θ in symmetric axis (Case 1)
The flow details of the liquid droplet and air for Case 1 are depicted in Fig.6，where the left side shows the dimensionless temperature field while the right side the flow structure.When the droplet interacts with the cold surface initially，the vortex is produced at the rim of the spreading front.With the evolution of the droplet morphology，the initial vortex gradually departs the rim of the droplet and an opposite vortex is produced due to the surface tension of the droplet，as shown in Fig.6（b）.When this vortex evolves with the liquid-air interface，some secondary vortexes are produced near the droplet interface att*=2.5.When the droplet enters the oscillation stage，the vortices merges，and the developed vortex，which is highlighted in Fig.6（d），gradually dissipates and disappears in the air，as shown in Fig.6（f）.
Fig.6 Evolutions of the temperature field (left side) and flow details (right side) for Case 1
Besides this，the temperature filed evolves from the bottom and shows the tendency to close the droplet，which is also indicated in Fig.5.
3.2 Effect of impact velocity
The droplets with different impact velocity collide on the clean aircraft wings will cause different iceshapes.In this section，the freezing of the droplet with velocity of 5 m/s to 20 m/s are simulated and the evolutions of the droplet morphologies and ice growing processes inside the droplet are shown in Fig.7.
When the impact velocity is small（Fig.7（a）），the droplet spreads less due to the low inertia and freezing of the bottom.However，the freezing is not fast enough to restrict the droplet lift from the surface when the surface tension of the droplet transits the potential energy to the kinetic energy as shown att*=1.5 andt*=2.5.When the droplet spreads on the ice after its top gets to the highest position，the second spreading causes the second ice layer att*=4.0.
With the increase of the impact velocity（Fig.7（b）），the inertia promotes the spreading and thus，the heat transfer between the droplet and the cold surface is enhanced.The momentum is rapidly lost compared with that for Case 5 so that the droplet forms an ice pancake att*=4.0 until it is totally frozen.For the droplet with impact velocity of 20 m/s（Case 7），the freezing is further enhanced.The first ice rim is formed att*=0.8 while the droplet is still in the spreading stage and the unfrozen liquid part of the droplet starts to form the second ice rim when the droplet spreads to its extreme att*=1.5，which is highlighted in Fig.7（c）.
Fig.7 Evolutions of droplet morphologies and ice growing processes for Cases 5 to 7
The heat transfer at the surface determines the ice growth.The temperature distribution along the radial direction at the bottom of the droplet is analyzed as presented in Fig.8.The temperature difference between the droplet bottom and the cold surface indicates the gratitude of the temperature which is related to the thermal flux density and reflects the strength of the heat transfer.It is noted that the largest temperature difference is near the spreading front of the droplet.For the droplet with small impact velocity（Case 5），the temperature of the droplet bottom keeps decreasing and that of the air near the surface keeps decreasing trend too.
Fig.8 Distribution of the dimensionless temperature θ in radial direction for Case 5
The evolutions of the temperature and the flow details are shown in Fig.9.For small impact velocity（Case 5），the vortex is produced at the rim of the droplet.With the droplet retraction and rebounding，the vortex develops and its range becomes the largest at the moment when the droplet rebounds to the highest position.And the secondary vortex is produced correspondingly which gradually evolves att*=2.5 as shown in the highlighted red box in Fig.9（a）.For Case 6，due to the rapid freezing of the droplet，the vortex produced at the droplet rim dissipates with time att*=4.0 as shown in Fig.9（b）.It is noted that a pair of vortices are produced at the droplet rim for Case 7.This is attributed to the freezing of the droplet rim，which is not smooth so that the flow structure is complicated here.Another finding is that with the increase of the impact velocity，the temperature of the air far from the droplet rim evolves slower.This could be explained by the characteristic time for heat transfer and dynamics.Despite in the same dynamic dimensionless time，the time of the heat transfer is still short for the droplet with larger impact velocity.
Fig.9 Evolutions of the temperature field (left side) and flow details (right side) for Cases 5 to 7
To further illustrate this，the relative magnitudes of the inertial force，viscous force and surface tension are evaluated byReandWenumbers.In Figs.10 and 11，the distributions of the dimensionless velocity in droplet symmetric axis defined asfor Case 3 and Case 7 are depicted，respectively.Reis 558.7 whileWevaries from 132.5 to 264.9.It’s obvious that the magnitudes of the maximum velocity in the symmetric axis are decreasing for both Cases 3 and 7.This is due to the viscous dissipation as well as the freezing process，which inhibits the dynamics of the air.
Fig.10 Distribution of the dimensionless velocity U*z in symmetric axis for Case 3 (Re=558.7,We=132.5)
Fig.11 Distribution of the dimensionless velocity U*z in symmetric axis for Case 7 (Re=558.7,We=264.9)
During the droplet impact and freezing on the cold surface，the momentum and heat transfer between the droplet and air is complicated.To analyze the droplet-air interactions，the dimensionless temperature distributions in the vertical direction atr=2Dfor Cases 5 and 6 are shown in Fig.12.It’s found that with the increase ofRefrom 139.7（Case 5）to 419（Case 6），the temperature becomes higher.For Case 6，the droplet spreads faster in the radial direction because of the larger inertial force，as shown in Fig.7（b）so that the release of the latent heat from the droplet during the freezing process weakens the cooling of the air from the cold surface.Therefore，the temperature is higher compared with that in Case 5.
Fig.12 Distribution of the dimensionless temperature in vertical direction at r=2D for Cases 5 and 6 (t*=0.5 and t*=4.0).
3.3 Effect of cold surface
The conditions for the aircraft icing include the micro-sized supercooled droplets，the low ambient temperature and so on.This section studies the effect of the cold surface，the temperature of which ranges from -10 ℃ to -25 ℃.The simulation results of the droplet morphological evolutions for Cases 8 to 10 are shown in Fig.13.
In dimensionless time，the dynamic behavior shows the similar characteristic.The droplet spreads and retracts to rebound.The frozen root of the droplet hinders the detachment from the cold surface.When the temperature of the cold surface is not low enough（Case 8），the droplet retracts and forms an intruded tip att*=4.0 as shown in Fig.13（a）.When the surface is cold enough（Case 10），the rapid freezing limits the dynamic behavior of the droplet，the center of the droplet is not able to restore to the flat shape and a slight crater is observed att*=4.0 as indicated in Fig.13（c）.
Fig.13 Evolutions of droplet morphologies and ice growing processes for Cases 8 to 10
To analyze the effect of the cooling from the cold surface，the temperature distributions along the droplet symmetric axis att*=4.0 for Cases 8 to 10 are shown in Fig.14.With the temperature decreasing on the cold surface，the dimensionless temperature inside the droplet at the same moment drops more quickly，indicating that the heat transfer inside the droplet is not proportional to the temperature of the cold surface.The increased temperature difference enhances the heat transfer so that the dimensionless temperature inside the droplet becomes lower.
Fig.14 Distribution of the temperature in symmetric axis at t* =4.0 for Cases 8 to 10
The micro-size droplet impact and freezing on the cold surface is simulated by the developed numerical method.The flow field is solved by the multiphase lattice Boltzmann flux solver，and the interface of the liquid and air is captured by the phase field method while the liquid-ice interface is tracked by the heat transfer of the three phases where the finite difference method is applied.The freezing characteristics of the droplet are investigated and the effects of the droplet size，impact velocity and the temperature of the cold surface are considered.The main conclusions are drawn as follows.
（1）The developed method is validated by the comparison between the computational result and the experimental data.It turns out that the coupled numerical method is able to accurately capture the dynamic behavior and freezing process of the droplet.
（2）For small droplets，the freezing from the droplet bottom is not able to completely inhibit the dynamic process，and the top shows the tendency to detach from the cold surface.While for lager-sized droplets，the freezing is enhanced and it shows a hatlike shape.The larger impact velocity causes the increased contact area，and the heat transfer is rapid，which restricts the droplet detachment from surface.The temperature of the cold surface governs the dynamic behavior of the droplet center.When the surface is cold enough，the droplet shows a crater in the center.
（3）The temperature inside the droplet in the symmetric axis shows a nonlinear distribution in the ice region.For the cold surface with different temperatures，the colder the surface is，the lower the dimensionless temperature becomes and the increased temperature difference accounts for the enhanced freezing.Therefore，the heat flux near the liquid-ice interface is smaller for the colder surface.
