水滴撞击冻结过程的数值模拟研究

边庆勇; 朱程香; 赵宁; 朱春玲; BIAN Qingyong; ZHU Chengxiang; ZHAO Ning; ZHU Chunling

Transactions of Nanjing University of Aeronautics & Astronautics

+Advanced Search

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Numerical Investigation on Droplet Impact and Freezing PDF

- ORCID：
BIAN Qingyong ¹
- ORCID：
ZHU Chengxiang ^1,2
- ORCID：
ZHAO Ning ^1,2
- ORCID：
ZHU Chunling ^1,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

CLC： V211.3

Updated：2023-05-08

DOI：10.16356/j.1005-1120.2023.02.007

OUTLINE

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.

Keywords

aircraft icing; supercooled droplet; impact; freezing; numerical simulation

0 Introduction

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^{［Reference 3-4}3-4］. Besides small droplets， larger-sized 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^{［Reference 5

Baidu Scholar}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^{［Reference 6

Baidu Scholar}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.^{［Reference 11

Baidu Scholar}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.^{［Reference 13

Baidu Scholar}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.^{［Reference 14

Baidu Scholar}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.^{［Reference 15

Baidu Scholar}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.^{［Reference 16

Baidu Scholar}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.^{［Reference 17

Baidu Scholar}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.^{［Reference 21-23}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.^{［Reference 24

Baidu Scholar}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.^{［Reference 25-27}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.^{［Reference 20

Baidu Scholar}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.^{［Reference 28-29}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.^{［Reference 30

Baidu Scholar}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.

1 Numerical Methods

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 liquid-ice 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

\frac{\partial p}{\partial t} + ρ c_{s}^{2} \nabla \cdot u = ρ c_{s}^{2} C (1 - \frac{ρ_{l}}{ρ_{s}}) \frac{\partial γ}{\partial t}

(1)

\begin{matrix} \frac{\partial ρ u}{\partial t} + \nabla \cdot (ρ u u) = \\ - \nabla p + \nabla \cdot [μ (\nabla u + {(\nabla u)}^{T})] + F_{s} \end{matrix}

(2)

where t is the time， p the pressure， ρ the density， u the velocity， c_s the speed of sound， μ the dynamic viscosity， and F_s the surface tension. The subscripts “l” and “s” represent the variables of liquid and ice， respectively. C is 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

\frac{\partial C}{\partial t} + \nabla \cdot (u C) = C (1 - \frac{ρ_{l}}{ρ_{s}}) \frac{\partial γ}{\partial t}

(3)

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 tension F_s is related to C， and the following form of the surface tension is applied^［

31-32］

F_{s} = - C \nabla μ_{C}

(4)

where μ_C is the chemical potential and can be calculated as^［

33］

μ_{C} = 2 β C (C - 1) (2 C - 1) - k_{f} \nabla^{2} C

(5)

where β and k_f are related to the surface tension σ and the width of the interface W in the simulation， yielding β = 12σ/W and k_f = 3σW/2. For the wettability of the cold surface， the static contact angle model is applied， which is described as^［

31， 34］

n \cdot \nabla C = - c o s φ^{e q} \frac{4 C_{w} (1 - C_{w})}{W}

(6)

where n is unit vector of the surface， C_w the non-air fraction at the cold surface and φ^eq the 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］

\frac{\partial ρ H}{\partial t} + \nabla \cdot (ρ u H) = \nabla \cdot (k \nabla T)

(7)

where H is the enthalpy of the fluid， T the temperature and k the thermal conductivity. During the freezing process of the droplet， the latent heat L is released and it is incorporated in the enthalpy H， which is defined by^［

35］

H = C_{p} T + C L γ

(8)

where C_p is the specific heat capacity. Due to the difference of the multiphase material properties， ρ， C_p and k should transit smoothly from liquid to air， ice to liquid and ice to air. Therefore， the material properties can be taken by

ψ = ψ_{a} + C [ψ_{s} + γ (ψ_{l} - ψ_{s}) - ψ_{a}]

(9)

where ψ represents ρ， C_p or k. 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 temperature T and expressed as follows^［

10， 35］

γ = \frac{1}{2} t a n h [k_{T} (T + T_{s h i f t})] + \frac{1}{2}

(10)

where k_T is the parameter to control the temperature range of phase transition and increase the numerical stability of the computation. Here， T_shift is used to ensure that the liquid temperature is at freezing point and matches the corresponding initial liquid fraction γ₀ = 1-C_p∆T/L， where ∆T is 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 fraction C and liquid fraction γ， the different phases could be identified. 1-C， Cγ and C（1-γ） denote air， liquid and ice， respectively^［

37］. In the ice region， the velocity is corrected by

u = u [1 - C (1 - γ)]

(11)

Obviously， when C（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^{［Reference 39-40}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^{［Reference 40

Baidu Scholar}40］.

2 Model Validation

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 of U and the cold surface is kept with constant temperature T_w while 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 diameter D is 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^{[Reference 41

Baidu Scholar}41]

Phase	Liquid	Air	Ice
ρ/(kg∙m^-3)	1 000	1.293	917
C_p/(J∙kg^-1∙K^-1)	4 216	1 005	2 060
k/(W∙m^-1∙K^-1)	0.569	0.024 4	2.2
μ/(Pa∙s)	1.79×10^-3	1.72×10^-5	N/A
L/(J∙kg^-1)	333 624	N/A	N/A
σ/(N∙m^-1)	0.075 5	N/A	N/A

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 number Re = ρ_lUD/μ_l and the Weber number We=ρ_lU²D/σ for this condition are Re=1 110.6， We=18.4， respectively， which are within the range of droplet impact condition in aircraft icing^{［Reference 1

Baidu Scholar}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^{［Reference 42

Baidu Scholar}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 diameter D to the initial droplet diameter D₀. 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]

3 Results and Discussion

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 ∆T is equal to the absolute value of T_w， i.e.， ∆T= |T_w|. 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

Case No.	Diameter D/μm	Velocity U/(m⋅s^-1)	Surface temperature T_w/℃
1	20	10	-15
2	50	10	-15
3	100	10	-15
4	200	10	-15
5	50	5	-15
6	50	15	-15
7	50	20	-15
8	50	10	-10
9	50	10	-20
10	50	10	-25

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 time t* and dimensionless temperature θ， which are defined as

t * = \frac{U t}{D}

(12)

θ = - \frac{T}{T_{w}}

(13)

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 10¹ μm to 10² μ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.

Fig.4 Evolutions of morphologies of different-sized droplets and ice growing processes for Cases 1 to 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 at t* = 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 at t* = 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 at t* = 4.0 as shown in Fig.4（c）.

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. At t* = 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 at t*=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 at t*=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.

Fig.7 Evolutions of droplet morphologies and ice growing processes for Cases 5 to 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 at t*=1.5 and t*=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 at t*=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 at t*=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 at t*=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 at t*=1.5， which is highlighted in Fig.7（c）.

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 at t*=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 at t*=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 by Re and We numbers. In Figs.10 and 11， the distributions of the dimensionless velocity in droplet symmetric axis defined as $U_{z}^{*}$ = U_z/U₀ for Case 3 and Case 7 are depicted， respectively. Re is 558.7 while We varies 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 at r=2D for Cases 5 and 6 are shown in Fig.12. It’s found that with the increase of Re from 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.

Fig.13 Evolutions of droplet morphologies and ice growing processes for Cases 8 to 10

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 at t* = 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 at t* = 4.0 as indicated in Fig.13 （c）.

To analyze the effect of the cooling from the cold surface， the temperature distributions along the droplet symmetric axis at t* = 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

4 Conclusions

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 hat-like 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.

Contributions Statement

Mr. BIAN Qingyong developed the numerical method， simulated the dynamic freezing process of the droplet， analyzed and discussed the results， and wrote the original draft. Dr. ZHU Chengxiang gave useful suggestions on the analysis of the results and reviewed the original draft. Prof. ZHAO Ning contributed to the freezing model and reviewed the original draft. Prof. ZHU Chunling reviewed and edited the original draft， and contributed to the discussion and background of the study. All authors commented on the manuscript draft and approved the submission.

Acknowledgements

This work was supported by the National Natural Science Foundation of China （No.11832012） and the National Major Scientific Research Instrument Development Project （No.12227802）.

Conflict of Interest

The authors declare no competing interests.

References

CAO Y， WU Z， SU Y， et al. Aircraft flight characteristics in icing conditions‍［J］. Progress in Aerospace Sciences， 2015， 74： 62-80. [Baidu Scholar]

ZHANG C，LIU H.Effect of drop size on the impact thermodynamics for supercooled large droplet in aircraft icing［J］.Physics of Fluids，2016，28（6）： 062107. [Baidu Scholar]

CAO Y， XIN M. Numerical simulation of supercooled large droplet icing phenomenon： A review‍［J］. Archives of Computational Methods in Engineering， 2019， 27（4）： 1231-1265. [Baidu Scholar]

DING L， CHANG S. Evaluation of inflight deicing performance for designed electrothermal systems with experimental verification‍［J］. Journal of Aircraft， 2020， 57（1）： 156-166. [Baidu Scholar]

YI X， WANG Q， CHAI C， et al. Prediction model of aircraft icing based on deep neural network‍［J］. Transactions of Nanjing University of Aeronautics and Astronautics， 2021， 38（4）： 535-544. [Baidu Scholar]

WANG Y， HAN L， ZHU C， et al. Design of an experimental set-up concerning interfacial stress to promote measurement accuracy of adhesive shear strength between ice and substrate‍［J］. Transactions of Nanjing University of Aeronautics and Astronautics， 2022， 39（5）： 561-568. [Baidu Scholar]

CHANG S， QI H， ZHOU S， et al. Experimental and numerical study on freezing process of water droplets under surfaces with different wettability‍［J］. Applied Thermal Engineering， 2023， 219： 119516. [Baidu Scholar]

MA Z， XIONG W， CHENG P. 3D lattice Boltzmann simulations for water droplet’s impact and transition from central-pointy icing pattern to central-concave icing pattern on supercooled surfaces. Part Ⅱ： Rough surfaces‍［J］. International Journal of Heat and Mass Transfer， 2021， 172： 121153. [Baidu Scholar]

RAKOTONDRANDISA A， DANAILA I， DANAILA L. Numerical modelling of a melting-solidification cycle of a phase-change material with complete or partial melting‍［J］. International Journal of Heat and Fluid Flow， 2019， 76： 57-71. [Baidu Scholar]

ATTARZADEH R， DOLATABADI A. Icephobic performance of superhydrophobic coatings： A numerical analysis‍［J］. International Journal of Heat and Mass Transfer， 2019， 136： 1327-1337. [Baidu Scholar]

LI W， WANG J， ZHU C， et al. Numerical investigation of droplet impact on a solid superhydrophobic surface‍［J］. Physics of Fluids， 2021， 33（6）： 063310. [Baidu Scholar]

ASHOKE RAMAN K， JAIMAN R K， LEE T S， et al. Dynamics of simultaneously impinging drops on a dry surface： Role of impact velocity and air inertia‍［J］. Journal of Colloid and Interface Science， 2017， 486： 265-276. [Baidu Scholar]

ZHANG B， SANJAY V， SHI S， et al. Impact forces of water drops falling on superhydrophobic surfaces‍［J］. Physics Review Letters， 2022， 129（10）： 104501. [Baidu Scholar]

ZHOU W， LONEY D， FEDOROV A G， et al. Lattice Boltzmann simulations of multiple-droplet interaction dynamics‍［J］. Physics Review E， 2014， 89（3）： 033311. [Baidu Scholar]

PALACIOS J， HERNáNDEZ J， GóMEZ P， et al. Experimental study of splashing patterns and the splashing/deposition threshold in drop impacts onto dry smooth solid surfaces‍［J］. Experimental Thermal and Fluid Science， 2013， 44： 571-582. [Baidu Scholar]

EBRAHIM M， ORTEGA A. Identification of the impact regimes of a liquid droplet propelled by a gas stream impinging onto a dry surface at moderate to high Weber number‍［J］. Experimental Thermal and Fluid Science， 2017， 80： 168-180. [Baidu Scholar]

GUO J， ZOU S， LIN S， et al. Oblique droplet impact on superhydrophobic surfaces： Jets and bubbles‍［J］. Physics of Fluids， 2020， 32（12）： 122112. [Baidu Scholar]

TEMBELY M， DOLATABADI A. A comprehensive model for predicting droplet freezing features on a cold substrate‍［J］. Journal of Fluid Mechanics， 2018， 859： 566-585. [Baidu Scholar]

BLAKE J， THOMPSON D， RAPS D， et al. Simulating the freezing of supercooled water droplets impacting a cooled substrate‍［J］. AIAA Journal， 2015， 53（7）： 1725-1739. [Baidu Scholar]

CHAUDHARY G， LI R. Freezing of water droplets on solid surfaces： An experimental and numerical study‍［J］. Experimental Thermal and Fluid Science， 2014， 57： 86-93. [Baidu Scholar]

YAO Y， LI C， TAO Z， et al. Experimental and numerical study on the impact and freezing process of a water droplet on a cold surface‍［J］. Applied Thermal Engineering， 2018， 137： 83-92. [Baidu Scholar]

YAO Y， LI C， ZHANG H， et al. Modelling the impact， spreading and freezing of a water droplet on horizontal and inclined superhydrophobic cooled surfaces‍［J］. Applied Surface Science， 2017， 419： 52-62. [Baidu Scholar]

YAO Y， YANG R， LI C， et al. Investigation of the freezing process of water droplets based on average and local initial ice fraction‍［J］. Experimental Heat Transfer， 2019， 33（3）： 197-209. [Baidu Scholar]

WANG Y， JU L， HAN D， et al. Numerical investigation of the impacting and freezing process of a single supercooled water droplet‍［J］. Physics of Fluids， 2021， 33（4）： 042114. [Baidu Scholar]

VU T V. Numerical study of solidification of a drop with a growth angle difference‍［J］. International Journal of Heat and Fluid Flow， 2020， 84： 108599. [Baidu Scholar]

VU T V， NGUYEN C， KHANH D. Direct numerical study of a molten metal drop solidifying on a cold plate with different wettability‍［J］. Metals， 2018， 8（1）： 47. [Baidu Scholar]

VU T V， TRYGGVASON G， HOMMA S， et al. Numerical investigations of drop solidification on a cold plate in the presence of volume change‍［J］. International Journal of Multiphase Flow， 2015， 76： 73-85. [Baidu Scholar]

ZHANG X， LIU X， MIN J， et al. Shape variation and unique tip formation of a sessile water droplet during freezing‍［J］. Applied Thermal Engineering， 2019， 147： 927-934. [Baidu Scholar]

ZHANG X， LIU X， WU X， et al. Simulation and experiment on supercooled sessile water droplet freezing with special attention to supercooling and volume expansion effects‍［J］. International Journal of Heat and Mass Transfer， 2018， 127： 975-985. [Baidu Scholar]

STITI M， CASTANET G， LABERGUE A， et al. Icing of a droplet deposited onto a subcooled surface‍［J］. International Journal of Heat and Mass Transfer， 2020， 159： 120116. [Baidu Scholar]

WANG Y， SHU C， HUANG H B， et al. Multiphase lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio‍［J］. Journal of Computational Physics， 2015， 280： 404-423. [Baidu Scholar]

WANG Y， SHU C， YANG L M. An improved multiphase lattice Boltzmann flux solver for three-dimensional flows with large density ratio and high Reynolds number‍［J］. Journal of Computational Physics， 2015， 302： 41-58. [Baidu Scholar]

LEE T， LIU L. Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces‍［J］. Journal of Computational Physics， 2010， 229（20）： 8045-8063. [Baidu Scholar]

YANG L， SHU C， YU Y， et al. A mass-conserved fractional step axisymmetric lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio‍［J］. Physics of Fluids， 2020， 32： 103308. [Baidu Scholar]

TEMBELY M， ATTARZADEH R， DOLATABADI A. On the numerical modeling of supercooled micro-droplet impact and freezing on superhydrophobic surfaces‍［J］. International Journal of Heat and Mass Transfer， 2018， 127： 193-202. [Baidu Scholar]

CAO Y， TAN W， WU Z. Aircraft icing： An ongoing threat to aviation safety‍［J］. Aerospace Science and Technology， 2018， 75： 353-385. [Baidu Scholar]

LYU S， WANG K， ZHANG Z， et al. A hybrid VOF-IBM method for the simulation of freezing liquid films and freezing drops‍［J］. Journal of Computational Physics， 2021， 432： 110160. [Baidu Scholar]

BIAN Q， ZHU C， WANG J， et al. Numerical investigation on the characteristics of single droplet deformation in the airflow at different temperatures‍［J］. Physics of Fluids， 2022， 34（7）： 073307. [Baidu Scholar]

LIU X D， OSHER S， CHAN T. Weighted essentially non-oscillatory schemes‍［J］. Journal of Computational Physics， 1994， 115（1）： 200-212. [Baidu Scholar]

SHU C W， OSHER S. Efficient implementation of essentially non-oscillatory shock-capturing schemes‍［J］. Journal of Computational Physics， 1988， 77（2）： 439-471. [Baidu Scholar]

LIU X， MIN J， ZHANG X， et al. Supercooled water droplet impacting-freezing behaviors on cold superhydrophobic spheres‍［J］. International Journal of Multiphase Flow， 2021， 141： 103675. [Baidu Scholar]

ZHANG X， LIU X， WU X， et al. Impacting-freezing dynamics of a supercooled water droplet on a cold surface： Rebound and adhesion‍［J］. International Journal of Heat and Mass Transfer， 2020， 158： 119997. [Baidu Scholar]

SHINAN C， LIANG D， MENGJIE S， et al. Numerical investigation on impingement dynamics and freezing performance of micrometer-sized water droplet on dry flat surface in supercooled environment‍［J］. International Journal of Multiphase Flow， 2019， 118： 150-164. [Baidu Scholar]