Investigation on Phase Transition and Collection Characteristics of Non-Spherical Ice Crystals with Eulerian and Lagrangian Methods

Abstract

Ice crystal icing occurs in jet engine compressors, which can severely degrade jet engine performance. In this paper, two different numerical calculation methods, the Eulerian method and the Lagrangian method, were used to evaluate the dynamics, mass transfer, heat transfer, phase transition and trajectory of ice crystals. Then, we studied the effects of initial diameter, initial sphericity, initial temperature of ice crystal, and relative humidity of airflow on the phase transition and collection characteristics of ice crystal particles. Results indicate that the non-spherical characteristics of ice crystals have a significant impact on their impingement limits and collection characteristics. The collection coefficient of unmelted ice crystals is positively correlated with the initial particle diameter and sphericity, and negatively correlated with the initial particle temperature and the relative humidity of airflow. The melting rate of ice crystal particles on the impact surface increases exponentially with the initial diameter of the particles, linearly increases with the relative humidity of the airflow and initial temperature of the particles, and exponentially decreases with the sphericity of the particles.

1. Introduction

Icing phenomenon including supercooled droplet icing and ice crystal icing is always a threat to flight safety in aviation. Supercooled droplet icing refers to the process in which spherical supercooled liquid droplets impact to the cold surface of an aircraft to form ice accumulation. And ice crystal icing is aimed at solid ice particles which can have all kinds of shapes and with a warm environment. Previous studies have always focused on the icing of supercooled droplets. However, in the last decade a growing number of research efforts have been undertaken to understand ice crystal icing. In 2006, Mason et al. [1] analysed the causes of numbers of aero-engine failures and found that most of these accidents were caused by ice crystal icing rather than supercooled droplets icing. The occurrence of ice crystal icing phenomenon in the core of the engine can lead to a reduction in thrust, engine surge and other serious effects. Furthermore, it can cause fatal accidents from the shedding of ice blocks. Research has shown that the presence of liquid water is necessary for ice crystals icing. The liquid water can be either supercooled droplets contained in clouds or produced by melting ice crystals. In the front section of the engine’s compressor there is a warmer region with a total temperature higher than the melting temperature [2]. Ice crystals will be partially melted after following the airflow into that region. And the melted water creates the conditions for the ice crystals’ adhesion, thus the ice crystal icing phenomenon occurs. Research by the National Research Council of Canada (NRC) has shown [3,4,5,6] that the most severe ice crystal icing occurs when the melting rate is 10~25%. There are various irregular shapes of ice crystal particles in the atmosphere. Different shapes have different drag coefficients and convective heat transfer coefficients, which will affect the trajectory of ice crystals and the phase transition process. The melting state and collection of ice crystals are extremely difficult to measure experimentally inside the engine. Thus, it is necessary to investigate the phase change characteristic and collection characteristic of ice crystals during their movement in warm environments numerically.

The methods for numerical investigation of phase transition characteristic and collection characteristic of ice crystals are mainly divided into Eulerian and Lagrangian methods. Ayan et al. [7] investigated the motion and phase transition of ice crystal particles based on the Lagrangian method in 2D. Grift et al. [8] studied the ice crystal trajectories of a turbofan compressor based on the Lagrangian method, but they had not considered the non-spherical characteristics of the ice crystals. Aouizerate et al. [9] used Lagrangian method to simulate flight accretion under mixed-phase conditions. Zhengzhi W [10] et al. numerically simulated ice crystal trajectories based on Lagrange method. But they did not consider the non-spherical characteristics and phase transition characteristics of ice crystals. Koichiro H [11] et al. numerically simulated the trajectories and impacts of ice crystals of an axial fan based on Lagrange method. Iuliano et al. [12] analysed the effect of the shape of the ice crystals on their trajectory and the heat transfer performance based on Eulerian method. Villedieu and Trontin et al. [13,14] calculated ice crystal impact properties and ice crystal icing under mixed-phase conditions based on Eulerian method. Norde et al. [15] calculated the trajectories of ice crystals in an engine based on the Eulerian and Lagrangian methods, respectively, but they neglected the density change of ice crystal particles during the phase transition.

In this paper, the phase transition and collection characteristics of ice crystals have been investigated in warm environments. Under the premise of considering the non-spherical properties of ice crystals, a motion model considering the phase transition of ice crystals have been established, and a method for calculating the local collection coefficients of the impact surface under the phase transition conditions have been developed, based on the Eulerian and Lagrangian methods, respectively. Meanwhile, the motion phase transition process of the NACA0012 airfoil and the local collection coefficient on the surface have been calculated. In addition, the results obtained by two different numerical methods are compared and discussed in detail. Finally, the effects of initial particle size, initial temperature, initial sphericity of particles and relative humidity of air flow on the characteristics of ice crystal phase transformation and local collection coefficient of impact surface have been analyzed.

2. Ice Crystal Models

2.1. Drag Coefficient Model

Particle motion is affected by a variety of forces such as drag, gravity, buoyancy, virtual mass, basset history force, and lift. For traditional supercooled droplets generally only needs to consider its drag. But, since the size of ice crystals is much larger than traditional supercooled droplets, this paper considers both drag and gravity, in contrast, the size of other forces can be ignored.

The non-spherical characteristics of ice crystals determine that their drag coefficient model is different from that of supercooled droplets. There are many scholars who have studied the drag coefficient of non-spherical particles through experiments and obtained the drag coefficient relation equation by data fitting method. Typical relational equations commonly used include: “Haider and Levenspiel Model [16]”, “Ganser Model [17]”, “Holzer and Sommerfeld Model [18]”, and “Xianzhi-Song Model [19]” (

Table 1. Ice crystal drag coefficient model.

Model	Expression
Haider and Levenspiel	$\begin{array}{l}{C}_d=\frac{24}{Re}\left(1+\mathrm{A}R{e}^B\right)+\frac{C}{1+DR{e}^{-1}},Re<2.5\times {10}^4\\ with,\\ A=\exp \left(2.3288-6.4581\ast \Phi +2.4486\ast {\Phi }^2\right)\\ B=0.0964+0.5565\ast \Phi \\ C=\exp \left(4.905-13.8944\ast \Phi +18.4222\ast {\Phi }^2-10.2599\ast {\Phi }^3\right)\\ D=\exp \left(1.4681+12.2584\ast \Phi -20.7322\ast {\Phi }^2+15.8855\ast {\Phi }^3\right)\end{array}$
Ganser	$\begin{array}{l}{C}_d=\frac{24}{Re{K}_1}\left(1+0.1118{\left(Re{K}_1{K}_2\right)}^{0.6567}\right)+\frac{0.4305K_2}{1+\frac{3305}{Re{K}_1{K}_2}}\\ with,\\ Isometric \text{:}\\ K_1={\left(\frac{1}{3}+\frac{2}{3}{\Phi }^{-1/2}\right)}^{-1}-2.25\frac{d}{d_{tube}},K_2={10}^{1.8148{\left(-\log \Phi \right)}^{0.5743}}\\ Non-Isometric:\\ K_1={\left(\frac{1}{3}{{\Phi }_{\perp }}^{-0.5}+\frac{2}{3}{\Phi }^{-1/2}\right)}^{-1}-2.25\frac{d}{d_{tube}},K_2={10}^{1.8148{\left(-\log \Phi \right)}^{0.5743}}\end{array}$
Holzer and Sommerfeld	$C_d=\frac{8}{Re}\frac{1}{\sqrt{{\Phi }_{\parallel }}}+\frac{16}{Re}\frac{1}{\sqrt{\Phi }}+\frac{3}{\sqrt{Re}}\frac{1}{{\Phi }^{3/4}}+\frac{0.42}{{\Phi }_{\perp }}{10}^{0.4{(-{\log }_{10}\Phi )}^{0.2}}$
Xianzhi-Song	$C_d=\frac{24}{Re\left({\Phi }^{0.65}{\Phi }_{\perp }^{0.3}\right)}\ast {\left(1+0.35\ast Re\right)}^{0.44}$
Clift and Gauvin	$C_d=\left[\frac{24}{Re}\ast \left(1+0.15\ast R{e}^{0.687}\right)\right]+\frac{0.42}{1+\left(42500/R{e}^{1.16}\right)}$

In the above formulas, $\Phi$, ${\Phi }_{\perp }$ and ${\Phi }_{\parallel }$ represent the sphericity, crosswise sphericity, and lengthwise sphericity, respectively.$Re$ refers to relative Reynolds number. $K_1$ and $K_2$ are defined as the Stokes’ shape factor and the Newton’s shape factor, respectively. $d$ and $d_{tube}$ refer to the particle equivalent diameter and the diameter of the tube in which the particle is settling, respectively.

). In this paper, the local collection coefficients of NACA0012 airfoils under different drag models are calculated by Eulerian and Lagrangian methods, respectively, according to the calculation conditions of the benchmark test case (

Table 2. Benchmark tese case calculation conditions.

Parameter	Unit	Value
$\mathrm{Pressure} (p_a$)	$\left[\mathrm{P}\mathrm{a}\right]$	94,900
$\mathrm{Temperature} (T_a$)	$\left[\mathrm{K}\right]$	293
$\mathrm{Velocity} (u_a$)	$\left[\mathrm{m}/\mathrm{s}\right]$	67
$\mathrm{Relative} \mathrm{humidity} (RH$)	$[\%]$	70
$\mathrm{Angle} \mathrm{of} \mathrm{attack} (AOA$)	$\left[\text{°}\right]$	2
$\mathrm{Ice} \mathrm{crystal} \mathrm{temperature} (T$)	$\left[\mathrm{K}\right]$	270
$\mathrm{Spherical} (\Phi$)	$\left[\text{-}\right]$	0.4

). It is compared with the widely used relation for the drag coefficient of spherical particles proposed by Clift and Gauvin [20], and the results are shown in Figure 1.

As can be seen from the Figure 1, the horizontal axis represents the ratio of the arc length from the local position to the leading edge of the airfoil to the chord length of the airfoil, $s$/$c$. And the positive values of $s$/$c$ represent the upper surface of airfoil, the negative values represent the lower surface, and values equal to 0 represent the leading edge of the airfoil. The results obtained by the two numerical calculation methods are basically consistent. Moreover, the impingement limits and local collection coefficients of “Holzer and Sommerfeld Model”, “Ganser Model”, “Haider and Levenspiel Model” are basically the same. However, the results of “Xianzhi-Song Model” are close to those of “Clift and Gauvin Model”, but different from those of other three models. The difference of the local ice crystal collection coefficient at the leading edge is about 5.99%, and the difference of the upper impingement limit is about 59.24%, the lower impingement limit is about 27.66%. Due to the lack of experimental data, it is impossible to judge the accuracy of these drag coefficient models, and in the subsequent work of this paper, the “Ganser Model” is chosen because it is relatively simple and the most widely used.

2.2. Convective Heat Transfer Model

The non-spherical characteristics of ice crystals also determine the convective heat transfer characteristics different from supercooled droplets. However, the effect of particle shape on convective heat transfer is much less studied than the effect of particle shape on the drag coefficient. Due to the difficulty of experimental measurement, numerical and analog methods have been used to study the convective heat transfer characteristics of non-spherical particles. Typical Nusselt number relations include: “Richter and Nikrityuk Model [21]”, ”Comer and Kleinstreuer Model [22]”and “Villedieu Model [13]” (

Table 3. Ice crystal convective heat transfer model.

Model	Expression
Richter and Nikrityuk	$Nu=1.76+0.55\Phi P{r}^{1/3}R{e}^{1/2}{\Phi }_{\perp }^{0.075}+0.0014P{r}^{1/3}R{e}^{2/3}{\left(\frac{\Phi }{{\Phi }_{\perp }}\right)}^{7.2}$
Comer and Kleinstreuer	$Nu=1.393R{e}^{0.348}exp\left(0.248\left(1-E\right)\right)$
Villedieu	$Nu=2\ast {\Phi }^{1/2}+0.55\ast P{r}^{1/3}\ast {\Phi }^{1/4}\ast R{e}^{1/2}$
Ranz and Marshall	$Nu=2+0.6\ast P{r}^{1/3}\ast R{e}^{1/2}, for Re\in \left[0,250\right]$

In the above formulas, $Nu$ refers to Nusselt number, $Pr$ refers to Prandtl number, and $E$ is defined as aspect ratio.

). In this paper, the local collection coefficients of NACA0012 airfoils under different convective heat transfer models are calculated by Eulerian and Lagrangian methods, respectively, according to the calculation conditions of the benchmark test case (Table 2). It is compared with the widely used relation for the convective heat transfer model of spherical particles proposed by Ranz and Marshall [23], and the results are shown in Figure 2.

As can be seen from the Figure 2, the results obtained by the two numerical calculation methods are basically consistent. Moreover, the impingement limits and local collection coefficients of “Richter and Nikrityuk Model”, “Comer and Kleinstreuer Model” and “Villedieu Model” are basically the same. However, the results of the spherical “Ranz and Marshall Model” differ significantly from those of the other three models. The main differences are as follows: the local collection coefficient of unmelted ice crystals in the area near the leading edge and the lower wing surface decreases to a certain extent, and the local collection coefficient of melted water increases to a certain extent, and the error at the leading edge is about 10.04%. This indicates that the Nusselt number obtained by the “Ranz and Marshall Model” is larger than the other three non-spherical models, resulting in a high melting rate, which is consistent with the conclusion of Noder [24]. In the subsequent work of this paper, the “Villedieu Model” is chosen because it is the most widely used.

2.3. Phase Transition Model

The temperature of the ice crystal particles will change during the movement due to the warm airflow around, and the energy exchange with the airlow is calculated by the following formula:

$$Q_{air}=\frac{\pi \ast d_p\ast Nu\ast k_{air}}{\Phi }\ast \left(T_{air}-T_p\right),$$

(1)

where, $k_{air}$ represents the heat conductivity coefficient of air, $d_p$ refer to the particle equivalent diameter, $T_{air}$ and $T_p$ refer to the temperature of air and particle, respectively.

According to the temperature of ice crystals, the phase transition process of ice crystals is divided into three stages:

The first stage, when the particle temperature is below the melting temperature, the particles are pure solid ice crystals, which will undergo sublimation or deposition, and the temperature continues to rise until it reaches the melting temperature.

The second stage, when the particle temperature is equal to the melting temperature, the particles are a mixture of ice and water, which will undergo melt and evaporation or condensation. It is also believed that the phase transition phenomena all start from the outermost side, and that the melting water film surrounds the inner ice crystal nucleus and will not shed from it.

The third stage, when the particle temperature is higher the melting temperature, the particles are pure liquid droplets, which will undergo evaporation or condensation. Due to the continuous heat exchange with the surrounding air, the temperature of the particles continues to rise until they reach dynamic equilibrium or completely evaporate.

The shape evolution of non-spherical ice crystal during phase transition process is shown in Figure 3. In the first stage, there is sublimation or deposition on the surface of the solid ice crystal, and in the second and third stages, there is evaporation or condensation on the surface of the water film. In this paper, the processes of sublimation, deposition, evaporation and condensation are described by mass diffusion equation. The mass diffusion phenomenon is caused by the difference of vapor content between the ice crystal surface and the air flow environment. According to the evaporation model proposed by Crowe et al. [25], the mass flux applied to icing conditions with low evaporation rate is as follows:

$${\dot{m}}_{ev}=\pi \ast {\rho }_a\ast D_v\ast d_p\ast Sh\ast \left(Y_{v,s}-Y_{v,\infty }\right),$$

(2)

where, $S_h$ is Sherwood number, its expression is:

$$Sh=2\ast {\Phi }^{1/2}+0.55\ast R{e}^{1/2}S{c}^{1/3}{\Phi }^{1/4},$$

(3)

where, $S_c$ is Schmidt number, its expression is:

$$Sc=\frac{{\mu }_a}{{\rho }_a{D}_v},$$

(4)

where, $D_v$ is the vapor diffusivity in the air, the empirical formula proposed by Schirmer is adopted:

$$D_v=D_0{\left(\frac{T_a}{T_0}\right)}^{1.81}\ast \frac{p_0}{p_a},$$

(5)

where, $D_0=2.26\ast {10}^{-5}{m}^2/s, T_0=273.15K, p_0={10}^{5}Pa$, and $p_a$ refers to air pressure.

$Y_{v,s}$ and $Y_{v,\infty }$ represent the mass fraction of the particle surface and the steam in the free flow, respectively. And the steam mass fraction is expressed as follows:

$$Y_v=\frac{m_v}{m_a+m_v}\approx \frac{m_v}{m_a}=\frac{M_v{X}_v}{M_a{X}_a}=\frac{M_v{p}_v}{M_a{p}_a},$$

(6)

where, $m_v$ and $m_a$ represent the quality of water vapor and air in the air stream, respectively. $M_v$ and $M_a$ represent molar mass of water vapor and air, respectively. $X_v$ and $X_a$ represent the volume fraction of water vapor and air, respectively. $p_v$ refers to vapor pressure.

When we calculate $Y_{v,s}$, $p_v$ refers to the vapor pressure $p_{v,s}$ on the surface of the particle, it can be considered as the saturated vapor pressure corresponding to the particle temperature $T$, $p_{v,s}=p_{sat}\left(T\right)$. According to the saturation vapor pressure formula proposed by Sonntag [26], it can be obtained:

$$p_{sat,i}=exp\left(-6024.53T^{-1}+29.327+1.06\ast {10}^{-2}T-1.32\ast {10}^{-5}{T}^2-0.4938lnT\right),$$

(7)

$$p_{sat,w}=exp\left(-6096.94\ast T^{-1}+21.24-2.71\ast {10}^{-2}T+1.674\ast {10}^{-5}{T}^2+2.4335lnT\right),$$

(8)

where, $p_{sat,i}$ and $p_{sat,w}$ refer to the saturated vapor pressure on the ice surface and water film surface, respectively.

When we calculate $Y_{v,\infty }$,$p_v$ refers to the vapor pressure of free flow $p_{v,\infty }$, it can be considered as: $p_{v,\infty }=R_h\ast p_{sat}\left(T_{air}\right)$. Where, $R_h$ refers to relative humidity. $p_{sat}\left(T_{air}\right)$ is the saturated vapor pressure corresponding to the ambient temperature $T_{air}$ of the air flow around the ice crystal particles.

In the second stage, the ice crystals will melt, and the melt water film will evaporate or condense. In that, the heat absorbed by the ice crystals from the environment is entirely used for the latent heat of melting of the ice crystals and the latent heat of evaporation or condensation of the water film. According to the energy balance, the following relationship is obtained:

$$Q_{air}={\dot{m}}_{ev}{L}_{ev}+{\dot{m}}_{melt}{L}_m,$$

(9)

where, ${\dot{m}}_{ev}$ and ${\dot{m}}_{melt}$ represent evaporation flux and melt flux, respectively. $L_{ev}$ and $L_m$ represent the latent heat of evaporation and melting, respectively.

With the mass diffusion equation, we can derive the melting mass flux of ice crystals:

$${\dot{m}}_{melt}=\frac{1}{L_m}\left(\pi \ast d_p\ast Nu\ast k_{air}\ast \left(T_{air}-T_m\right)/\Phi -{\dot{m}}_{ev}{L}_{ev}\right),$$

(10)

In the first stage of the phase transition process, the particles always remain pure solid ice crystals, the sphericity maintains the initial sphericity, and the density is that of solid ice. In the third stage, the particles always remain pure liquid droplets with sphericity of 1 and density of liquid water. In the second phase, the particle is a mixed solid-liquid state, and the calculation of sphericity and density depends on the solid-liquid ratio of the mixed state, the calculation formula is as follows:

$$\begin{array}{l}Eulerian:\rho =\frac{{\alpha }_i{\rho }_i + {\alpha }_w{\rho }_w}{{\alpha }_i + {\alpha }_w}\\ Eulerian:\Phi =\left(\frac{{\alpha }_i}{{\alpha }_i + {\alpha }_w}\right)\ast {\Phi }_{ini} + \left(\frac{{\alpha }_w}{{\alpha }_i + {\alpha }_w}\right)\ast 1\end{array},$$

(11)

$$\begin{array}{l}Lagrangian:\rho =\frac{m_{p\_i}+m_{p\_w}}{m_{p\_i}/{\rho }_i + m_{p\_w}/{\rho }_w}\\ Lagrangian:\Phi =\left(\frac{m_{p\_i}/{\rho }_i}{m_{p\_i}/{\rho }_i + m_{p\_w}/{\rho }_w}\right)\ast {\Phi }_{ini} + \left(\frac{m_{p\_w}/{\rho }_w}{m_{p\_i}/{\rho }_i + m_{p\_w}/{\rho }_w}\right)\end{array},$$

(12)

where, ${\Phi }_{ini}$ refers to the initial sphericity of ice crystals. ${\rho }_i$ and ${\rho }_w$ refer to the density of ice and water, respectively.

3. Numerical Simulation

3.1. The Assumptions

In this paper, the following assumptions are proposed in the calculation of ice crystal icing problem: (1) Ice crystal particles don’t rotate during their motion and are always in the most stable dynamic state; (2) Ice crystal particles don’t collide or fuse during their motion; (3) The density of unmelted ice crystals and melted water remains constant; (4) During the melting process, particles remain at a constant melting temperature.

3.2. Eulerian Method

Eulerian method is calculated by solving a set of partial differential equilibrium equations based on the conservation law, which considers that particles are continuous phases, so the derivation of Eulerian equation involves the mean value in the phase space. In terms of computational cost, Eulerian method is generally more efficient. Its cost does not depend on the number of particles present in the computational configuration particularly. However, dealing precisely with complex physical phenomena, such as the crossing of two particles carrying jets or the interaction between particles and walls, Eulerian method is more cumbersome.

The phase transition characteristics of ice crystal particles will affect the mass conservation equation, momentum conservation equation and energy conservation equation of ice crystal particles at the same time. So, the source term is added to the governing equation to handle it in this paper. Since the melting of ice crystals will produce liquid water, this paper describes the change of volume fraction of melted water by adding a new equation. The complete eulerian control equations for ice crystal particles in three stages of phase transition are shown in

Table 4. Eulerian governing equations for non-spherical ice crystals.

Stage	Item	Equation
1	Mass	$\frac{\partial {\alpha }_i{\rho }_i}{\partial t}+\nabla ·\left({\alpha }_i{\rho }_i\overrightarrow{u}\right)=-{\dot{\alpha }}_{sub}{\rho }_i=-\frac{6{\rho }_a{D}_{v}Sh}{\Phi {\rho }_i{d}_p^2}\left(Y_{v,s}-Y_{v,\infty }\right)\alpha {\rho }_i$
	Mass	$\frac{\partial {\alpha }_w{\rho }_w}{\partial t}+\nabla ·\left({\alpha }_w{\rho }_w\overrightarrow{u}\right)=0$
	Momentum	$\frac{\partial \alpha {\rho }_p\overrightarrow{u}}{\partial t}+\nabla ·\left(\alpha {\rho }_p\overrightarrow{u}\overrightarrow{u}\right)=\frac{3{\mu }_a{C}_{d}Re}{4d_p^2}\alpha \left(\overrightarrow{u_a}-\overrightarrow{u}\right)-{\dot{\alpha }}_{sub}{\rho }_p\overrightarrow{u}$
	Energy	$\frac{\partial \alpha {\rho }_{p}T}{\partial t}+\nabla ·\left(\alpha {\rho }_{p}T\overrightarrow{u}\right)=\frac{6k_{a}Nu}{\Phi c_i{d}_p^2}\alpha \left(T_a-T\right)-{\dot{\alpha }}_{sub}{\rho }_i\left(\frac{L_{sub}}{c_i}+T\right)$
2	Mass	$\frac{\partial {\alpha }_i{\rho }_i}{\partial t}+\nabla ·\left({\alpha }_i{\rho }_i\overrightarrow{u}\right)=-\frac{6}{\Phi {\rho }_p{d}_p^2{L}_{melt}}\left(k_{a}Nu\left(T_a-T\right)-{\rho }_a{D}_v{L}_{ev}Sh\left(Y_{v,s}-Y_{v,\infty }\right)\right)$
	Mass	$\frac{\partial {\alpha }_w{\rho }_w}{\partial t}+\nabla ·\left({\alpha }_w{\rho }_w\overrightarrow{u}\right)=\frac{{\rho }_i}{{\rho }_w}\ast {\dot{\alpha }}_w\ast {\rho }_w-\frac{6{\rho }_a{D}_{v}Sh}{\Phi {\rho }_w{d}_p^2}\ast \left(Y_{v,s}-Y_{v,\infty }\right)\ast \alpha {\rho }_w$
	Momentum	$\frac{\partial \alpha {\rho }_p\overrightarrow{u}}{\partial t}+\nabla ·\left(\alpha {\rho }_p\overrightarrow{u}\overrightarrow{u}\right)=\frac{3{\mu }_a{C}_{d}Re}{4d_p^2}\alpha \left(\overrightarrow{u_a}-\overrightarrow{u}\right)-{\dot{\alpha }}_{ev}{\rho }_p\overrightarrow{u}$
	Energy	$\frac{\partial \alpha {\rho }_{p}T}{\partial t}+\nabla ·\left(\alpha {\rho }_{p}T\overrightarrow{u}\right)=\left(\frac{{\rho }_i}{{\rho }_w}-1\right){\dot{\alpha }}_{melt}{\rho }_p{T}_{melt}-{\dot{\alpha }}_{ev}{\rho }_p\left(\frac{L_{ev}}{c_w}+T_{melt}\right)$
3	Mass	$\frac{\partial {\alpha }_i{\rho }_i}{\partial t}+\nabla ·\left({\alpha }_i{\rho }_i\overrightarrow{u}\right)=0$
	Mass	$\frac{\partial {\alpha }_w{\rho }_w}{\partial t}+\nabla ·\left({\alpha }_w{\rho }_w\overrightarrow{u}\right)=-{\dot{\alpha }}_{ev}{\rho }_w=-\frac{6{\rho }_a{D}_{v}Sh}{\Phi {\rho }_w{d}_p^2}\left(Y_{v,s}-Y_{v,\infty }\right)\alpha {\rho }_w$
	Momentum	$\frac{\partial \alpha {\rho }_p\overrightarrow{u}}{\partial t}+\nabla ·\left(\alpha {\rho }_p\overrightarrow{u}\overrightarrow{u}\right)=\frac{3{\mu }_a{C}_{d}Re}{4d_p^2}\alpha \left(\overrightarrow{u_a}-\overrightarrow{u}\right)-{\dot{\alpha }}_{ev}{\rho }_p\overrightarrow{u}$
	Energy	$\frac{\partial \alpha {\rho }_{p}T}{\partial t}+\nabla ·\left(\alpha {\rho }_{p}T\overrightarrow{u}\right)=\frac{6k_{a}Nu}{\Phi c_w{d}_p^2}\alpha \left(T_a-T\right)-{\dot{\alpha }}_{ev}{\rho }_w\left(\frac{L_{ev}}{c_w}+T\right)$

In the above formulas, $c_i$ and $c_w$ represent specific heat capacity at constant pressure of ice and water, respectively. ${\alpha }_i$, ${\alpha }_w$ and $\alpha$ represent volume fraction of unmelted ice crystals, melted water and total particles, respectively. ${\dot{\alpha }}_{sub}$, ${\dot{\alpha }}_{melt}$ and ${\dot{\alpha }}_{ev}$ represent the rates of sublimation, melting and evaporation of volume fractions, respectively.

.

3.3. Lagrangian Method

The calculation of Lagrange method is to track the movement of particles in the flow field, establish the particle motion equation on the basis of the calculation results of the air flow field, solve the particle trajectory, and obtain the impact characteristics through a large number of particle trajectory calculation. Therefore, its accuracy directly depends on the number of sampled particles, and in general, convergence is slow. According to the central limit theorem, the convergence order is $n_p^{-1/2}$, where $n_p$ is the number of sampledparticles. As the number of particles increases, the computational cost doubles. Mean-while, the Lagrangian method is slightly more complex in determining the particle impingement limit, as well as the number and spacing of released particles. In this paper, a Lagrange method for calculating the Impingement limit trajectory Is proposed as follows (taking the two-dimensional case as an example):

(1)

Firstly, find two trajectories, one bypassing the upper surface of the frozen target and the other bypassing the lower surface of the frozen target. The horizontal coordinate of the starting point of the trajectory corresponds to the horizontal coordinate of the entrance of the calculation domain, and the vertical coordinate of the starting point is denoted as $y_u$ and $y_d$, respectively;

(2)

Determine whether the particle with a vertical coordinate of $y_m=\left(y_u+y_d\right)/2$ at the entrance of the computational domain collides with the icing target. If the collision occurs, then the starting ordinate of the upper limit trajectory is between $y_u$ and $y_m$, and the starting ordinate of the lower limit trajectory is between $y_m$ and $y_d$;

(3)

Further determine whether the trajectory with ordinate $y=\left(y_u+y_m\right)/2$ at the entrance of the calculation domain collides with the icing target. If it collides, then $y_m=y$; if the trajectory bypasses the upper surface, then $y_d=y$;

(4)

Repeat step (3) until the value of $y$ reaches convergence, and the upper limit particle trajectory can be obtained. And the trajectory of the lower limit particle can be obtained using the same method;

(5)

If the particle with a vertical coordinate of $y_m=\left(y_u+y_d\right)/2$ in step (2) does not collide with the frozen target, then based on whether it bypasses the object from top or bottom, use the methods described in steps (3) and (4) to find a new trajectory that intersects with the frozen target, with its starting point longitudinally marked as $y_m$, and then execute steps (3) and (4).

Based on the Lagrangian method, the motion equation, mass equation, and temperature equation of particles are obtained through the analysis of particle forces, mass transfer, and heat transfer, respectively. Based on convective heat transfer model and phase transition model, the mass transfer and heat transfer terms caused by particle phase transition are added to the equation to make the equations closed and balanced. The complete Lagrange governing equations of the three stages in the phase transition process of ice crystal particles are shown in

Table 5. Lagrangian governing equations for non-spherical ice crystals.

Stage	Item	Equation
1~3	Motion	$m_p\frac{d{\overrightarrow{u}}_p}{dt}=\frac{1}{8}{C}_D{\rho }_a\pi d_p^2\left({\overrightarrow{u}}_a-{\overrightarrow{u}}_p\right)\left|{\overrightarrow{u}}_a-{\overrightarrow{u}}_p\right|+\frac{1}{6}\pi d_p^3\left({\rho }_p-{\rho }_a\right)\overrightarrow{g}$
1	Mass	$\frac{d{m}_p}{dt}=-{\dot{m}}_{sub}=-\pi {\rho }_a{D}_v{d}_{p}Sh\left(Y_{v,s}-Y_{v,\infty }\right)$
	Temperature	$m_p{c}_i\frac{d{T}_p}{dt}=\pi d_p\frac{Nu}{\Phi }{k}_a\left(T_a-T_p\right)-{\dot{m}}_{sub}{L}_{sub}$
2	Mass	$\frac{d{m}_i}{dt}=-{\dot{m}}_{melt}=-\frac{1}{L_{melt}}\left(\pi d_{p}Nu{k}_{air}\left(T_{air}-T_{melt}\right)/\Phi -{\dot{m}}_{ev}{L}_{ev}\right)$
	Mass	$\frac{d{m}_w}{dt}={\dot{m}}_{melt}-{\dot{m}}_{ev}$
	Momentum	$T_p=T_{melt}$
3	Mass	$\frac{d{m}_p}{dt}=-{\dot{m}}_{ev}=-\pi {\rho }_a{D}_v{d}_{p}Sh\left(Y_{v,s}-Y_{v,\infty }\right)$
	Mass	$m_p{c}_w\frac{d{T}_p}{dt}=\pi d_p\frac{Nu}{\Phi }{k}_a\left(T_a-T_p\right)-{\dot{m}}_{ev}{L}_{ev}$

In the above formulas, $m_i$, $m_w$ and $m_p$ represent the mass of unmelted ice crystals, melted water and total particles, respectively. ${\dot{m}}_{sub}$, ${\dot{m}}_{melt}$ and ${\dot{m}}_{ev}$ represent the mass flux of sublimation, melting and evaporation, respectively.

.

3.4. Calculation of Local Collection Coefficient

In the icing problem, the local collection coefficient is the key parameter to calculate the impact water distribution along the surface, and it is a parameter to characterize the water collection capacity of a certain element surface. For supercooled droplets, in two-dimensional case, the local collection coefficient is calculated as follows: (three-dimension is similar)

$$Lagrangian:\beta =\frac{d{y}_{\infty }}{dy},$$

(13)

where, $d{y}_{\infty }$ represent the longitudinal distance of particles at the release position, and $dy$ represent the longitudinal distance of particles at the impact surface.

$$Eulerian:\beta =\frac{{\rho }_d\ast \alpha \ast \left(\overrightarrow{u}·\overrightarrow{n}\right)}{\left|u_{d\_\infty }\right|\ast LW{C}_{\infty }},$$

(14)

where, $\overrightarrow{u}$ refers to the velocity vector of particles. $\overrightarrow{n}$ refers to the normal vector of the surface of the local element. $u_{d\_\infty }$ refers to the particle velocity at the entrance. And $LW{C}_{\infty }$ refers to the liquid water content at the entrance.

Different from the supercooled droplets icing, ice crystals icing will undergo phase transition process during movement. The local total collection coefficient of ice crystals is divided into two parts, one is the collection coefficient of unmelted ice crystals ${\beta }_i$, and the other is the collection coefficient of melted water ${\beta }_w$. On the other hand, it is because of the existence of the phase transition process that the paper assumes that the unmelted ice crystals and melted water that hit the surface are all collected without bouncing and splashing. According to the definition, the formula for calculating the local collection coefficient is:

$$\begin{array}{l}Lagrangian:{\beta }_i=\frac{d{y}_{\infty } \ast IWC}{dy \ast IW{C}_{\infty }}\\ Lagrangian:{\beta }_w=\frac{d{y}_{\infty } \ast LWC}{dy \ast IW{C}_{\infty }}\end{array},$$

(15)

$$\begin{array}{l}Eulerian:{\beta }_i=\frac{{\rho }_i \ast {\alpha }_i \ast \left(\overrightarrow{u}·\overrightarrow{n}\right)}{\left|u_{d\_\infty }\right| \ast IW{C}_{\infty }}\\ Eulerian:{\beta }_w=\frac{{\rho }_w \ast {\alpha }_w \ast \left(\overrightarrow{u}·\overrightarrow{n}\right)}{\left|u_{d\_\infty }\right| \ast IW{C}_{\infty }}\end{array},$$

(16)

where, $IW{C}_{\infty }$ refers to the ice crystals content at the entrance. $IWC$ and $LWC$ are defined as unmelted ice crystals content and melted water content at the impact surface, respectively.

In the calculation of this paper, the ice crystal content is measured using dimensionless volume fraction. Taking the ice crystal content in the air at the entrance as the benchmark, the volume fraction of the ice crystal particles is dimensionless, that is, the volume fraction (dimensionless volume fraction) of the ice crystal particles at the entrance is 1, and the ice crystal content at other locations in the flow field is the relative volume fraction and the ratio of the ice crystal volume fraction at the entrance. If the particle content at any location in the flow field is required, it is only necessary to multiply the relative volume fraction of the particles at that location by the ice crystal content at the entrance of the flow field. The governing equations of the eulerian method directly contain the variable of the volume fraction, while the governing equations of the lagrangian method use the mass variable, so it is necessary to convert the mass variable of the ice crystal (or liquid water) into the volume fraction variable. Unless otherwise stated, the volume fractions mentioned below are all relative volume fractions.

4. Validation of Models and Methods

4.1. Calculation Conditions

Hauk [27] conducted 222 experiments on ice crystal particles suspended in an acoustic levator in order to study the phase transition process of ice crystal particles under forced convection environment. In this paper, considering the economic efficiency of calculations 27 phase transition processes with very clear experiment results are selected as the verification basis of which 12 are spherical particle and the remaining 15 are non-spherical particle. The selected test conditions include spherical and non-spherical particles, high relative humidity (36~74%) airflow and low relative humidity (2~4%) airflow. They all include different flow temperature, ice crystal initial diameter, initial temperature and initial sphericity for spherical and non-spherical particles. The flow rate, however, should not be too high (about 0.5~2 m/s) to ensure that the particle remains suspended in the acoustic field instead of being blown away. The detailed test conditions are shown in

Table 6. Test case on particle phase change: spherical particles.

Case	${\mathit{d}}_{\mathit{p}}$$\left[\mathbf{\mu }\mathbf{m}\right]$	$\mathbf{\Phi }\left[\text{-}\right]$	$\mathit{T}$$\left[\mathbf{K}\right]$	${\mathit{u}}_{\mathit{a}}$$\left[\mathbf{m}/\mathbf{s}\right]$	${\mathit{p}}_{\mathit{a}}$$\left[\mathbf{P}\mathbf{a}\right]$	${\mathit{T}}_{\mathit{a}}$$\left[\mathbf{K}\right]$	$\mathit{R}\mathit{H}[\%]$
1	715	1	256.4	1.00	94,900	293.2	4
2	994	1	256.4	1.00	94,900	293.2	4
3	915	1	257.0	1.25	97,220	293.2	4
4	1013	1	256.2	1.00	94,900	293.2	75
5	978	1	256.5	1.00	94,900	293.2	75
6	775	1	256.4	0.75	95,100	288.2	64
7	591	1	256.9	0.75	95,100	288.2	64
8	584	1	255.5	1.25	95,300	298.3	3
7	521	1	255.5	1.25	95,300	298.2	3
10	779	1	254.3	0.75	95,300	298.2	44
11	768	1	254.3	0.75	95,300	298.3	44
12	961	1	254.9	1.75	95,400	298.2	3

and

Table 7. Test case on particle phase change: non-spherical particles.

Case	${\mathit{d}}_{\mathit{p}}$$\left[\mathbf{\mu }\mathbf{m}\right]$	$\mathbf{\Phi }\left[\text{-}\right]$	$\mathit{T}$$\left[\mathbf{K}\right]$	${\mathit{u}}_{\mathit{a}}$$\left[\mathbf{m}/\mathbf{s}\right]$	${\mathit{p}}_{\mathit{a}}$$\left[\mathbf{P}\mathbf{a}\right]$	${\mathit{T}}_{\mathit{a}}$$\left[\mathbf{K}\right]$	$\mathit{R}\mathit{H}[\%]$
13	784	0.51	257.5	1.00	95,870	292.9	4
14	551	0.70	255.5	1.00	95,870	292.8	4
15	1071	0.84	257.0	1.00	94,900	293.1	4
16	690	0.49	256.0	0.75	95,300	288.3	61
17	1013	0.78	257.9	1.00	95,600	293.2	78
18	656	0.78	254.9	1.25	95,300	298.2	3
19	572	0.68	255.1	1.25	95,300	298.2	3
20	929	0.66	254.1	0.75	95,300	298.2	44
21	845	0.59	254.2	0.75	95,300	298.3	44
22	634	0.81	254.6	1.25	95,400	298.3	40
23	699	0.59	254.8	1.25	95,400	298.3	40
24	732	0.79	255.1	1.75	95,400	298.3	3
25	792	0.67	254.6	1.25	95,400	303.3	2.2
26	915	0.69	255.9	1.25	95,400	298.2	56
27	662	0.79	255.9	1.25	95,330	298.3	56

.

The verification domain is a two-dimensional square region with a size of 0.05 × 0.05 m. The boundary condition of the entrance on the left side of the calculation domain is the velocity entrance, the right side is the outflow boundary, and the upper and lower is the symmetric boundary. The airflow and ice crystal particles are one-way coupled, so the calculation of the airflow field is converged first. Then, the phase transition process of the ice crystal particles is calculated using the Eulerian method and the Lagrangian method, respectively.

4.2. Verification Results and Analysis

The phase transition process of ice crystals includes three stages. Figure 4 and Figure 5 show the changes in particle temperature over time for spherical particles (case 9 and 11) and non-spherical particles (case 25 and 26), respectively. Figure 6 and Figure 7 show the changes of volume fraction of unmelted ice crystal and melted water over time for spherical particles (case 9 and 11) and non-spherical particles (case 25 and 26), respectively. From the temperature perspective, particle temperature is divided into three stages, corresponding to the three stages of ice crystal transition. In the first stage, the ice crystal temperature always rapidly rises from the initial temperature to the melting temperature (273.15 K). The second stage mainly involves the melting of ice crystal particles, and the temperature is maintained at the melting temperature of 273.15 K. At this time, the ice crystal particles are partially melted and accompanied by evaporation or condensation of the water film. The third stage is the heating stage of the droplet. When the ice crystal particles are completely melted into the droplets, the droplets absorb heat from the warm environment and continues to warm up. Finally, the droplet reaches a steady-state temperature, and the heat absorbed by the droplet from the environment reaches a dynamic balance with the heat taken away by the evaporation effect.

From the perspective of volume fraction, there are also three stages. The volume fraction of solid ice crystal equal to 1, and the volume fraction of liquid water is 0, in the initial state. In the first stage, the ice crystal volume fraction of case 9 and case 25 showed an increasing trend, while the ice crystal volume fraction of case 11 and case 26 showed a decreasing trend, which was related to the relative humidity of the airflow and the temperature of the particles. In the first stage, the particle temperature is quite low. When the relative humidity of the airflow is also low, the vapor content on the surface of the ice crystals is greater than that in the airflow, leading to sublimation and a decrease in the volume fraction of the ice crystals; When the relative humidity of the airflow maintains higher, the vapor content on the surface of ice crystals is lower than that in the airflow, leading to condensation and an increase in the volume fraction of ice crystals. In addition, at this stage, the volume fraction of liquid water is always 0. However, the first stage of pure solid ice crystals quickly transitions to the second stage. In the second stage, ice crystal melting dominates, accompanied by water film evaporation or condensation phenomena. The volume fraction of unmelted ice crystals gradually decreases, while the volume fraction of melted water gradually increases. When ice crystals completely melt into droplets, the volume fraction of unmelted ice crystals becomes 0, and the volume fraction of melted water is greater than 0 but less than 1, the second stage ends. After reaching the third stage, the particles become pure liquid, and the volume fraction of melted water decreases over time, regardless of whether it is high or low relative humidity. This is because the temperature of the water droplet become higher, causing the vapor content on the surface of the water droplet to be greater than that in the airflow, resulting in evaporation rather than condensation. Moreover, the lower the relative humidity of the airflow, the faster the evaporation rate.

The non-spherical characteristics of particles mainly affect the resistance coefficient and heat transfer coefficient for ice crystal particles. For the trend of temperature changes and volume fraction changes during the three phase transition processes, the non-spherical characteristics of particles do not show any special features compared with spherical particles.

Ice crystal particle melting time and particle equivalent diameter at complete melting under the 27 conditions mentioned above have been calculated using Lagrangian and Eulerian methods respectively. And the calculated results have been compared with the experimental results of Hauk [27]. Figure 8 shows the comparison of ice crystal particle melting time between the calculated and experimental values including spherical particles and non-spherical particles. And Figure 9 shows the comparison of the particle equivalent diameter at complete melting. It can be seen from the figures that the calculated results are close to the experimental results under the conditions of different incoming temperature, incoming velocity, relative humidity, ice crystal initial diameter, initial temperature and initial sphericity. Figure 10 shows the errors between experimental and calculated values for melting time and equivalent diameter, respectively. The positive value of the vertical axis means that the calculated value is greater than the experimental value, and the negative value is opposite. The errors of melting time are mostly within 10% with the maximum errors are 19.20% for Lagrangian method and 18.09% for Eulerian method. And the errors of equivalent diameter are mostly within 2% with the maximum errors are 3.04% for Lagrangian method and 2.48% for Eulerian method.

The ice crystal particles produced in the experiment have a certain surface roughness, but numerical calculations did not consider this factor. However, this roughness may affect its surface heat transfer characteristics, which may be the main reason for the error in melting time between experiments and calculations. It may be necessary to improve the heat transfer models to take into account the effects of ice crystal surface roughness in the future. At the same time, there may be certain errors in measuring the sphericity and diameter of ice particles during the experimental process. The calculation results are roughly consistent with the experimental results. The physical model and calculation method used in this paper can effectively simulate the phase transition process of ice crystal particles in warm airflow environments. On the other hand, the calculation results of Eulerian method and Lagrangian method are basically consistent, which indicates that both methods are applicable to the problem studied in this paper. Also, the results of the two different numerical calculation methods can be mutually validated due to the little experimental data.

5. Calculation of Phase Transition Process and Local Collection Coefficient

5.1. Calculation Conditions

Apply the calculation model and method of ice crystal motion and phase transition established in the previous text to the calculation of the flow field around an airfoil, while calculating the local collection coefficient of the impact surface. The geometric model uses NACA0012 airfoil with a chord length of 1 m and a calculation domain of 10 m × 6 m × 6 m. After grid independence verification, the final total mesh quantity is about 80,000, and the bottom mesh height of the wing surface is 0.001 mm. The mesh near the wing is shown in Figure 11. The calculation conditions are as follows: the air flow temperature is 293 K, the pressure is 94,900 Pa, the velocity is 67 m/s, the relative humidity of the air flow is 70%, and Angle of attack is 2°. The average equivalent diameter of the ice crystal particles is 220 $\mu \mathrm{m}$, the initial temperature is 270 K, and the initial sphericity is 0.4.

5.2. Results and Discussion

The calculation results of the volume fraction of unmelted ice crystals and melted water near the wing are shown in Figure 12. In the results of two different numerical calculation methods, the obvious shadow zone can be seen in the trailing edge of the wing, the location of the impingement limit and the size of the shadow zone are very consistent, there is no ice crystal particles, and the volume fraction of solid ice crystal and liquid water are all equal to 0. Eulerian method calculates the volume fraction throughout the entire computational domain, so the region of increased unmelted ice crystal (melted water) concentration can be seen outside of the shadow zone. However, the Lagrangian method only releases particles within the impingement limit in order to reduce computational costs, and therefore cannot observe that region. As ice crystal particles move backwards, they continuously absorb heat from the external environment. And the melting process continues, resulting in a gradual decrease in the volume fraction of unmelted ice crystal and an increase in the volume fraction of melted water.

The local collection coefficients of unmelted ice crystals and melted water are shown in Figure 13. The local collection coefficients of the stagnation point are the largest, which are 0.75 and 0.11 respectively, and the melting rate there can be calculated to be 13.75%. The collection coefficient gradually decreases along the chord until the upper and lower impingement limits are reached. And the impingement limit of the lower wing surface is greater than that of the upper wing surface, because the angle of attack is 2°, which is consistent with the volume fraction contours. The results obtained by two different numerical methods, Eulerian method and Lagrangian method, are basically consistent.

5.3. Influence Parameter Analysis

The calculation conditions are shown in

Table 8. Calculation conditions of different influence parameters.

Run	${\mathit{d}}_{\mathit{p}}$$\left[\mathbf{u}\mathbf{m}\right]$	$\mathbf{\Phi }\left[\text{-}\right]$	$\mathit{T}$$\left[\mathbf{K}\right]$	$\mathit{R}\mathit{H}[\%]$
a	50	1.0	270	4
	70	1.0	270	4
	120	1.0	270	4
	170	1.0	270	4
b	170	0.4	270	70
	220	0.4	270	70
	320	0.4	270	70
	520	0.4	270	70
c	120	1.0	270	10
	120	1.0	270	40
	120	1.0	270	70
	120	1.0	270	100
d	120	1.0	240	70
	120	1.0	250	70
	120	1.0	260	70
	120	1.0	270	70
e	320	0.1	270	70
	320	0.4	270	70
	320	0.7	270	70
	320	1.0	270	70

. Calculation conditions of different influence parameters., the effects of initial diameter, initial sphericity, initial temperature of ice crystal, and relative humidity of airflow on the phase transition and collection characteristics of ice crystal particles has been analyzed, separately.

5.3.1. Effects of Ice Crystal Initial Equivalent Diameter

The influence of initial equivalent diameter of ice crystal on the local collection coefficient and the melting rate of the leading edge are shown in Figure 14, Figure 15 and Figure 16, respectively. As the initial equivalent diameter of ice crystal increases, the impingement limits of unmelted ice crystals and melted water increases. This is because the larger the diameter of ice crystals, the stronger the effect of ice crystals following the airflow, resulting in a greater impingement limit. Meanwhile, as the particle diameter increases, the collection coefficient of unmelted ice crystals gradually increases, while the collection coefficient of liquid water gradually decreases, and the melting rate of ice crystal particles generally decreases exponentially. This is because the heat required for ice crystal melting increases as the particle diameter increases. The heat transfer area on the particle surface also increases as the particle diameter increases, resulting in an increase in heat absorbed from the environment, but the rate of volume growth is greater than the rate of area growth. The larger the particle size, the smaller the relative surface area. Therefore, the amount of heat required due to volume increase is greater than the amount of heat absorbed due to area increase. On the other hand, the heat absorbed by ice crystal melting is latent heat, and the surface heat transfer of ice crystal particles is sensible heat. The latent heat of melting is much larger than the sensible heat of surface heat transfer.

5.3.2. Effects of Relative Humidity

The influence of relative humidity (RH) on the local collection coefficient and the melting rate of the leading edge are shown in Figure 17 and Figure 18, respectively. As the relative humidity of the flow increases, the impingement limits of unmelted ice crystals and melted water remain unchanged. The collection coefficient of unmelted ice crystals decreases gradually, while the collection coefficient of melted water increases gradually, and the melting rate of ice crystal particles increases linearly. In the second stage of the ice crystal transition process, the ice crystal will melt and evaporate or condense. From the perspective of energy balance, the latent heat of melting of ice crystal particles is equal to the sum of the heat absorbed by the particles from the airflow and the latent heat of evaporation or condensation. When the relative humidity is low, the water content in the airflow is small, and the liquid water tends to evaporate and absorb heat. Moreover, the higher the relative humidity is, the less the mass transfer by evaporation of liquid water will be. When the relative humidity is high, the water content in the flow is high, and the liquid water tends to condense and release heat. Therefore, as the relative humidity of the airflow increases, the latent heat of ice crystal melting increases, the melting rate increases, the collection coefficient of liquid water increases, and the collection coefficient of solid ice crystals decreases.

5.3.3. Effects of Ice Crystal Initial Temperature

The influence of ice crystal initial temperature on the local collection coefficient and the melting rate of the leading edge are shown in Figure 19 and Figure 20, respectively. As the initial temperature increases, the impingement limits of unmelted ice crystals and melted water remain unchanged. The collection coefficient of unmelted ice crystals decreases gradually, while the collection coefficient of melted water increases gradually, and the melting rate of ice crystal particles increases linearly. In the first stage of the ice crystal phase transition process, the ice crystal temperature rises from the initial temperature to the melting temperature of 273.15 K. The higher the initial temperature, the shorter the duration of the first stage, and the earlier the melting process begins. After transitioning to the second stage, the particle temperature is all 273.15 K, so the phase transition process of ice crystals is completely consistent in the second stage at different initial temperatures. Therefore, for the same movement distance, the earlier the melting process starts, the longer the melting process lasts, the greater the melting rate, the greater the collection coefficient of melted water, and the lower the collection coefficient of unmelted ice crystals.

5.3.4. Effects of Ice Crystal Initial Sphericity

The influence of ice crystal initial sphericity on the local collection coefficient and the melting rate of the leading edge are shown in Figure 21 and Figure 22, respectively. As the initial sphericity increases, the impingement limits of unmelted ice crystals and melted water increases. By analyzing the drag coefficient model of non-spherical ice crystals, it can be concluded that under the same other conditions, as the sphericity increases, the drag coefficient of ice crystals decreases. Therefore, the effect of ice crystals following the airflow become stronger, leading to a greater impingement limit.

At the same time, it can be seen from the figures that with the increase of sphericity, the collection coefficient of unmelted ice crystals increases, the collection coefficient of melted water decreases, and the melting rate of ice crystal particles generally decreases exponentially. The smaller the sphericity of the particle, the more irregular the shape of the particle. Therefore, the greater the surface area of particles of the same equivalent diameter, the efficiency of the particles in absorbing heat from the environment is enhanced. From the perspective of energy balance, the latent heat of ice crystal melting has been increased. So, the melting rate of ice crystals has increased, the local collection coefficient of unmelted ice crystals has decreased, and the local collection coefficient of melted water has increased.

6. Conclusions

In this article, we developed a computational model for the phase transition of ice crystal particles and a method for calculating the local collection coefficients of unmelted ice crystals and melted liquid water. A three-stage ice crystal phase transition control equation and calculation method were established based on the Eulerian and Lagrangian method, and the calculation model and methods were validated. At the same time, we compared and analyzed the two different numerical calculation methods. Finaly, based on NACA0012 airframe, the phase transition and collection characteristics of ice crystals were calculated, and the parameters affecting the collection characteristics of melted ice crystals were further calculated and analyzed. The following conclusions were obtained:

• In the benchmark case of this paper, the error of the leading edge collection coefficient between the non-spherical particle drag model selected in this paper and the widely used spherical particle drag model is 13.07%, the upper limit error is about 59.24%, and the lower limit error is about 27.66%. And the error of heat transfer model is about 10.04% at the leading edge of collection coefficient. Therefore, it is necessary to consider the non-spherical characteristics for ice crystal icing.
• The numerical calculation results of the phase transition model developed in the paper are basically consistent with the experimental results. The maximum error of melting time and particle equivalent diameter is 19.20% and 3.04%, respectively.
• The impingement limit and collection coefficient distribution obtained using two different numerical calculation methods are very consistent, without considering the rebound and splashing of ice crystals. Developing two different numerical calculation methods can mutually verify each other, due to the extreme lack of experimental data. Lagrangian method is relatively complex in determining the impingement limit and the number and spacing of released particles.
• The melting rate of ice crystal particles near the leading edge is mostly within the plateau period of high ice crystal icing severity (5~25%) for NACA0012. The results of ice crystal melting rate can be effectively applied to the risk assessment of ice crystal icing.
• Impingement limit of ice crystal particles increases with the increase of particle diameter, increases with the increase of sphericity, and does not change with the initial temperature or the relative humidity. The collection coefficient of unmelted ice crystals increases with the increase of particle diameter, decreases with the increase of relative humidity, decreases with the increase of initial particle temperature, and increases with the increase of initial particle sphericity. And the collection coefficient of melted water is the opposite.
• The melting rate of ice crystal particles increases exponentially with the particle diameter, linearly with the relative humidity of airflow, linearly with the initial temperature of the particles, and exponentially decreases with the sphericity of the particles.