Three-Dimensional Trajectory and Impingement Simulation of Ice Crystals Considering State Changes on the Rotor Blade of an Axial Fan

Abstract

Ice crystal icing occurs in jet engine compressors, which can severely degrade jet engine performance. In this study, we developed an ice crystal trajectory simulation, considering the state changes of ice crystals with a forced convection model, indicating a significant difference in impinging ice crystal content on the blade for tiny ice crystals. Then, ice crystal trajectory simulations were performed for the rotor blade of an axial fan to investigate the effects of ice crystal size and relative humidity on collision characteristics. The results indicate that the surrounding air affects the composition of tiny ice crystals before collision, and the flight time until impingement on the rotor blade varies significantly depending on the span position. Among them, ice crystals with a diameter of 50 μm impinge with water content that is most likely to adhere to the blade. Three-dimensional simulation results show that many ice crystals impinge not only on the leading edge, where icing occurs as revealed by the two-dimensional simulations but also on the trailing edge of the hub side. This study emphasizes the importance of evaluating the three-dimensional impingement position and water content in the prediction of ice crystal icing.

1. Introduction

Some clouds contain many super-cooled water droplets and ice crystals. When an aircraft flies through such clouds, the super-cooled water droplets impinge on the solid walls and form an ice layer. This icing phenomenon in a jet engine also poses a threat to safe navigation because it can deteriorate the aerodynamic performance by changing the fan blade and vane shapes. Furthermore, it can cause fatal accidents from the shedding of ice blocks. Based on the liquid water content in the atmosphere and ambient temperature, the ice accretion by super-cooled water droplets can occur in various kinds of icing, such as glaze and rime icing. For the prediction of ice shape, various icing models have been developed. The pioneering work in developing an icing model, which considers mass and energy conservation laws, was proposed by Messinger [1]. Myers [2,3] improved the Messinger model by considering ice formation accompanied by thin liquid film flow. Then, Özgen [4] developed the extended Messinger model by applying the two-dimensional extension of Myers’s model, expressed as a one-dimensional mathematical model. Other icing models have been proposed worldwide, such as the LEWICE code developed by the Glenn Research Center Center [5,6], the CAPTA code and IGLOO2D developed by ONERA [7,8], and TRAJICE by DRA [9]. In recent years, many icing models have been developed for super-cooled large droplet (SLD) icing [10,11,12] and for rime icing using the coupled simulation of grid-based and particle-based methods [13,14].

Icing simulation generally consists of four steps: grid generation, flow field simulation, droplet trajectory simulation, and thermodynamics computation. Long-term icing simulations were performed by repeating the computational procedure of these four steps (i.e., multi-shot simulation). Droplet trajectory simulation has been carried out using either Eulerian or Lagrangian methods. The Eulerian method is based on the continuous particle field transport equation and can apply the highly particle-laden flow, supported not only by one-way coupling [15,16,17], but also by two-way coupling [18]. Recently, the one-way coupling method has been incorporated into open CFD solvers [19], and Long et al. [20] applied the Eulerian wall film model for icing simulation using the Euler–Euler approach. In contrast, the Lagrangian method is based on the trajectory equation of individual particles, allowing for trajectory simulation with particle size distribution and individual complex wall impingement [12,21].

Meanwhile, ice crystals have been observed to bounce off walls and are generally not considered to cause icing even in sub-zero environments because they are solid particles. However, ice crystal icing has been reported in environments exceeding the freezing point [22]. This phenomenon is attributed to ice crystals existing in a low-pressure compressor with a mainstream static temperature of approximately 303 K (30 °C). Under these conditions, the icing factors and environment are significantly different from those of icing caused by super-cooled droplets. First, ice crystals melt as they fly through a high-temperature environment, which causes them to evaporate or condense into a water film. Next, the impingement of ice crystals on the wall surface involves complex phenomena, such as break-up, bounce-off, or sticking, depending on the impact speed, melting rate, and wetting, indicating that these factors must be considered when predicting ice crystal icing. In addition, it is difficult to observe the actual phenomena because ice crystal icing occurs inside the jet engine during flight. Therefore, the generating process of ice crystal icing is not yet fully understood.

Previous studies focused on some of the above-mentioned factors. Lou and Hammond [23] developed a heat and mass transfer model for ice crystals flying in a high-temperature environment. They estimated the temporal variation in the size of ice crystals sucked into the engine using this model. Palacios et al. [24] conducted experiments in which melted ice crystals impacted an aluminum wedge. The results showed that the ice accretion area on the surface expanded with increases in impact angle, impact velocity, and melting ratio. Yan and Palacios [25] floated a single ice crystal using an ultrasonic levitator and experimentally investigated the correlation between environmental temperature and the melting rate over time. Hauk [26] performed experiments on the impingement of ice crystals under various conditions. He considered four different break-up modes for impinging ice crystals on a dry aluminum flat plate and found that the velocity of the break-up threshold is higher for impingement on a steel sphere with a water film compared to the dry-wall condition (i.e., without a water film). Alvarez et al. [27] proposed an empirical model for determining the threshold velocity of ice crystal break-up by impinging stainless steel with ice crystals at various melting ratios. They also found that the threshold velocity required to break the ice crystals increased significantly as the melting ratio increased.

Experimental and numerical studies have been conducted on airfoils and their actual shapes to understand the phenomena in jet engines. Currie et al. [28] experimentally investigated the relationship between the ratio of liquid water content (LWC) to total water content (TWC) and the adhesion efficiency in mixed-phase icing. They found that at a Mach number of $0.25$, the sticking efficiency is maximum at an LWC/TWC of approximately 10–20%, indicating that icing tends to occur. Struk et al. [29] conducted mixed-phase icing experiments on a NACA0012 airfoil to investigate the effects of wet-bulb temperature, particle size distribution, and Mach number on icing shape. Aouizerate et al. [30] used CEDRE, developed by ONERA, to simulate ice crystal trajectories in an engine inlet and fan, and in a fan and low-pressure compressor. Their simulation aimed to estimate changes in particle size distribution due to fracturing and the proportions of water and ice in the flow field. Flegel [31] conducted experiments on the behavior of ice crystals passing through a jet engine flow path, investigating the change in ice crystal size distribution caused by break-up downstream of the engine fan. They also examined the effect of the upstream ice crystal size distribution on the surface temperature and accumulated mass at the stator vanes. The above phenomena (e.g., break-up, melting, and solidification of ice crystals) have been modeled by real-scale simulations, and the particle size and water content distributions of ice crystals have been analyzed. Nagorski et al. [32] performed Eulerian simulations of the flow field and particle trajectories in the first stage of NASA Stage 67 to investigate the boundary conditions leading to ice crystal icing. Their results highlighted the importance of fully reproducing the boundary layer, as the flow in this layer is crucial for adequately melting the ice crystals. Iwago et al. [33] performed ice crystal icing simulations on a two-dimensional compressor stator blade of a NACA65-210. Furthermore, they performed icing simulations, considering the temperature distribution on the surface and inside of the stator blade; the temperature inside the blade dropped below the freezing point at the leading edge where the ice crystals impinged, resulting in the formation of an ice layer. More recently, Neuteboom and Chalmers [34] developed experimental equipment to study ice crystal icing in aircraft engines. They observed the growth of ice crystal icing and found out that the melt ratio of the ice crystal is an important|parameter.

Thus, models for various factors related to the occurrence of ice crystal icing have been developed through experimental and numerical studies, and simulations have been conducted on actual machine shapes. However, these previous studies did not sufficiently investigate the three-dimensional mass and water content distribution of ice crystals impinging on the blades. The ice crystal composition affects ice crystal icing via the cooling rate inside the airfoil and water film formation; therefore, it is essential to clarify the ice crystal composition before collision. In addition, it is necessary to clarify the impingement properties of ice crystals on the rotor blade as well as the stator vane based on the compressor structure to improve the prediction accuracy of ice crystal icing. Therefore, in this study, we developed an ice crystal trajectory simulation in a three-dimensional flow field around a rotor blade. Then, we introduced a state change model, which includes melting, evaporation, and condensation of ice crystals in flight, to accurately reproduce actual icing phenomena. The simulations were performed for four cases of median volume diameter (MVD) and three cases of relative humidity, and the effect of applying the state change model on impingement characteristics was numerically evaluated. This study highlights the importance of considering all three dimensions when predicting ice crystal icing, as the impingement mass and water content vary greatly depending on the span and chord position, as well as ice crystal size and relative humidity conditions.

2. Ice Crystal Model

2.1. State Change Model of the Ice Crystal

Ice crystal particles vary in phase composition based on the heat received from the surrounding air during the flight. To compute the phase composition, we adopted the state change model [23], which considers the phase changes of melting, evaporation, and condensation of the ice crystal. The assumptions of the phase state model are as follows:

• The melted liquid phase completely and concentrically covers all spherical ice particles (i.e., one-dimensional approximation in the radial direction).
• Heat transfer inside an ice crystal occurs only via heat conduction (i.e., without convective flow in the water film), and heat transfer from the surrounding air to the ice crystal occurs via convective heat transfer while assuming a small ice crystal.
• Either evaporation or condensation occurs at the water–air interface, and latent heat is considered to accompany the phase change at the interface.

Figure 1 shows the schematics of the ice crystal model for a partially melted state during flight, composed of spherical ice of radius $r_i$ and a surrounding water film of external surface radius $r_w$; the interface temperature between ice and water film is consistent with the melting temperature of ice $T_{melt}$. We considered three heat transfers between the interfaces: $Q_1$, $Q_2$, and $Q_3$. Here, $Q_1$ represents the heat received from the surrounding air, $Q_2$ represents the heat transferred into the ice from the water film, and $Q_3$ represents the heat of evaporation or condensation [35] on the outer water film surface.

From the heat transfer estimation, the time evolution of the water film radius $r_w$, ice particle radius $r_i$, and water film surface temperature $T_w$ were evaluated by the ice crystal model based on that of Lou and Hammond [23], as follows:

$$\begin{array}{ccc}\hfill \frac{d{r}_w}{dt}& =& \frac{D_c{M}_w}{{\rho }_{w}R}\left(\frac{P_{s,\infty }}{T_{\infty }}RH-\frac{P_{s,w}}{T_w}\right)\frac{1}{r_w},\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \frac{d{r}_i}{dt}& =& -\left\{r_i{h}_m\left(T_{\infty }-T_w\right)+\frac{D_c{M}_w}{R}\left(\frac{P_{s,\infty }}{T_{\infty }}RH-\frac{P_{s,w}}{T_w}\right)L_{ev}\right\}\frac{r_w}{r_i^2{\rho }_i{L}_{melt}},\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill k_w{r}_i\frac{T_w-T_{melt}}{r_w-r_i}& =& r_i{h}_m\left(T_{\infty }-T_w\right)+\frac{D_c{M}_w}{R}\left(\frac{P_{s,\infty }}{T_{\infty }}RH-\frac{P_{s,w}}{T_w}\right)L_{ev},\hfill \end{array}$$

(3)

where $T_{\infty }$ and t represent the air temperature at the ice crystal location and time, respectively; R, $RH$, and $h_m$ represent the universal gas constant, relative humidity, and heat transfer coefficient at a film temperature of $T_m=(T_{\infty }+T_w)/2$, respectively, for the surrounding air; $M_w$, $k_w$, and ${\rho }_w$ denote the molecular mass of water, heat conductivity [36], and density of water, respectively; $L_{ev}$ and $L_{melt}$ represent the latent heat of evaporation of water [37] and melting of ice, respectively; and ${\rho }_i$ represents the density of ice. In addition, $P_{s,\infty }$ and $P_{s,w}$ represent the saturated vapor pressures at $T_{\infty }$ and $T_w$, respectively, and $D_c$ denotes the diffusion coefficient of the vapor in the surrounding air.

After the ice crystal was terminated to melt, i.e., became a droplet, the droplet diameter was estimated by Equation (1), and the temperature change in the droplet [38] was estimated as follows:

$$\begin{array}{ccc}\hfill W_w{C}_w\frac{d{T}_w}{dt}& =& 4\pi r_w^2{h}_m\left(T_{\infty }-T_w\right)+\frac{d{W}_w}{dt}{L}_{ev},\hfill \end{array}$$

(4)

where $C_w$ and $W_w$ represent the constant-pressure specific heat of water and mass of water, respectively. Additionally, the droplet temperature $T_w$ is defined as the mean value within the droplet.

The heat transfer coefficient $h_m$ in Equations (2)–(4) can be estimated by the heat convection model for a spherical particle. By assuming a sufficiently small particle, a constant Nusselt number model ($Nu=2$) can be adopted, as in the previous study by Lou and Hammond [23], where the Nusselt number is defined by the heat conductivity of air $k_a$ and particle diameter $2r_w$. In this study, we considered the Nusselt number dependency because our objective ice crystal was relatively large. Therefore, we evaluated $h_m$ using the following three models of the Nusselt number $Nu$: The Nusselt number without convection $N{u}_{WO}$, with natural convection $N{u}_N$ [39], and with forced convection $N{u}_F$ around a sphere are expressed as follows:

$$\begin{array}{ccc}\hfill N{u}_{WO}& =& 2,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill N{u}_N& =& 2+\frac{0.589R{a}^{\frac{1}{4}}}{{\left(1+{\left(\frac{0.469}{Pr}\right)}^{\frac{9}{16}}\right)}^{\frac{4}{9}}},\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill N{u}_F& =& 2+0.60P{r}^{\frac{1}{3}}R{e}_{ic\_m}^{\frac{1}{2}},\hfill \end{array}$$

(7)

respectively, where $Pr$ is the Prandtl number of surrounding air, and $R{e}_{ic\_m}$ is the particle Reynolds number based on the diameter of the ice crystal $d_{ic}(=2r_w)$, the viscosity of air ${\mu }_m$ is estimated by $T_m$, alongside the relative velocity between the surrounding flow and ice crystal. $Ra$ represents the Rayleigh number, which is defined as follows:

$$\begin{array}{ccc}\hfill Ra& =& \frac{g{\beta }_m(T_{\infty }-T_w)d_{ic}^3}{{\nu }_m{\alpha }_m},\hfill \end{array}$$

(8)

where g, ${\beta }_m$, ${\nu }_m$, and ${\alpha }_m$ denote acceleration due to gravity, volume coefficient of expansion, kinetic viscosity, and thermal diffusivity of water, respectively. In addition, the variables with subscript m represent the values estimated based on film temperature $T_m$. The constant parameters used in our simulation are summarized in

Table 1. Constant parameters for the state change model.

Melting point of water	$T_{melt}$	[K]	$273.15$
Density of water	${\rho }_w$	[kg/m${}^3$]	$999.0$
Density of ice	${\rho }_i$	[kg/m${}^3$]	$917.0$
Vapor diffusion coefficient	$D_c$	[m${}^2$/s]	$22.0\times {10}^{-6}$
Molecular mass of water	$M_w$	[kg/mol]	$18.02\times {10}^{-3}$
Universal gas constant	R	[J/(K·mol)]	$8.314$
Latent heat of melting of ice	$L_{melt}$	[kJ/kg]	$333.5$

.

2.2. Ice Crystal Melting

Convection affects the Nusselt number and melting of the ice crystal. Here, we compare the melting times of the ice crystal in the numerical simulation for three different models:

• Case A: State change model without convection (i.e., non-convection), $Nu$ is computed by Equation (5), and the water film temperature varies.
• Case B: State change model with natural convection, $Nu$ is computed by Equation (6), and the water film temperature varies.
• Case C: NASA model, the effect of natural convection is considered, but the water film temperature is uniform [40] and calculated according to Soltis [41].

For the computational conditions, we followed the experimental study by Soltis [41], where the ice crystals were floated by ultrasonic levitation technology, and the total melting time was measured from a frozen ice crystal with the melting temperature $T_{melt}$ to a completely melted state for 14 cases: ice crystal diameter $0.324$–$1.112$ mm, ambient temperature 289–304 K, and relative humidity 41–100%.

Figure 2 compares the melting times in the numerical simulations and from the experimental results [41]. The marker shows a single simulation, e.g., for the melting time of 15 s in the experiment, the numerical simulations predicted 25 s, 22 s, and 21 s for Cases A, B, and C, respectively. Therefore, the upper area of the diagonal line corresponds to the underestimation of the melting time and vice versa. The numerical simulation reasonably predicted the melting time, although there was a discrepancy in the experiment. At melting times from 6 s to 12 s in the simulation, the melting time was overestimated compared with the experimental results because the relative humidity during the experiment fluctuated significantly [41] and induced flow around the ice crystal by the acoustic levitation was not negligible for the small ice crystal melting [42,43]. For the longer melting time cases, Case C underestimated the melting time and Case A overestimated the melting time. In the former and latter, the effects of the temperature variation in the water film of the ice crystal and convection were not considered, respectively. In contrast, Case B reasonably and relatively reproduced the melting time because the model considers the natural convection and temperature distribution in the liquid film.

Accordingly, we found that both convection and temperature variations in the water film are important in the melting simulation of ice crystals. In the following, instead of natural convection, using Equation (6), we employed forced convection using Equation (7) to simulate the melting ice crystal process in an axial fan because the ice crystal flows through the passage.

3. Numerical Simulation

3.1. Target and Flow Conditions

We performed an ice crystal icing simulation for an axial fan (Showa Denki Co., Osaka, Japan, Kairyu series A2D6H-411), as shown in Figure 3. We considered a computational domain, including one blade, and the periodic boundary condition was imposed in the circumferential direction. An overset-grid system was employed, including a main grid and sub-grid, as shown in Figure 4. The grid point numbers on the main grid and the sub-grid were $151\times 62\times 101$ and $262\times 81\times 251$, respectively, and the total number of grid points was approximately 6,270,000. The validity of this grid system was confirmed in our previous studies [44,45]. The total temperature, total pressure, flow angle, and turbulent quantities were fixed as boundary conditions at the inlet, while the Mach number was extrapolated. The static pressure was specified at the outlet, and the other variables were extrapolated.

Table 2. Computational flow conditions [46].

Rotational speed	[rpm]	1800
Mass flow rate of air	[kg/s]	$5.21$
Inlet static temperature	[K]	$303.15$
Relative humidity ($RH$)	[%]	40, 60, 80
Ice water content (IWC)	[g/m${}^3$]	$7.0$
Median volume diameter (MVD)	[μm]	25, 50, 100, 200

summarizes the numerical conditions. The rotational speed and mass flow rate of air were set according to a previous experiment [46]. The other parameters were determined by referring to an environment with ice crystal icing [33]. In this study, we investigated the effect of a state change in ice crystals on rotor blade icing for different relative humidities (40%, 60%, and 80%) and MVDs (25, 50, 100, and 200 μm). One million ice crystals were injected at two chord lengths upstream of the rotor blade. The actual amount of impinging ice crystals is estimated by Parcel approximation based on the ice crystal condition of the MVD and ice water content (IWC) in Table 2. The initial temperature of the ice inside the ice crystal is set to be the melting temperature $T_{melt}$. The initial positions of the ice particles were randomly distributed on the cross-section of the flow passage, and the initial velocity was set as the local velocity. Regarding the number of injected ice crystal particles, it was confirmed that the impingement rate of injected ice crystals on the blade surface for the $1.0$ million particle injection case shows no difference compared to the $0.5$ and $0.75$ million particle injection cases.

3.2. Numerical Procedure

Numerical simulations of ice crystal trajectory considering a state change were performed based on the Euler–Lagrange method with one-way coupling. The airflow field was assumed to be frozen; the steady flow field was simulated. In addition, mesh reconstruction was not performed since we focused on the impinging ice crystal composition before the ice adhesion. The numerical procedure involved the following three steps: (1) generation of computational grids, (2) computation of the flow field, and (3) computation of the ice crystal trajectories and acquisition of the impingement characteristics of the ice crystals on the rotor blade surface. These steps are explained in the following.

The computational grid system employed an overset grid composed of the main grid for the entire flow path and a sub-grid around the blade, as shown in Figure 4. This grid system was generated using an elliptical hyperbolic differential equation.

The steady flow field was obtained by the Reynolds averaged numerical simulation (RANS) using the in-house code. The flow was assumed to be a three-dimensional compressible turbulent flow. The governing equations were the continuity, Navier–Stokes, and energy equations with the Kato–Launder k-ε [47] and high Reynolds number wall function models. Second-order upwind total variation diminishing (TVD) [48] and second-order central difference schemes were applied to the spatial discretization for the inviscid and viscous terms, respectively, based on the finite differential scheme. The LU-ADI scheme [49] was used for time integration.

The ice crystal trajectories were simulated by the Lagrangian approach. Each ice crystal was assumed to be a sphere without rotation, and break-up and collisions were neglected before impingement on the rotor blade. In addition, we assumed a low-volume fraction of ice crystals in the atmosphere, and we ignored the effect of turbulent fluctuation on the particle’s trajectory, such as clustering and transport, because the particle Stokes number is relatively large. Therefore, we used a one-way coupling method, i.e., ice crystals do not affect the flow field, but the flow affects the ice crystal motion. Hence, we employed the simplified Basset–Boussinesq–Oseen (B-B-O) equation as the governing equation for ice crystal trajectory computation, which assumes that only flow drag, centrifugal, and Coriolis forces act on the ice crystals, as follows:

$$\begin{array}{ccc}\hfill \frac{d{\overrightarrow{U}}_{ic}}{dt}& =& \frac{3}{4}{C}_D\frac{{\rho }_a}{{\rho }_{ic}}\frac{1}{d_{ic}}{\overrightarrow{U}}_r\left|{\overrightarrow{U}}_r\right|-\left\{2\overrightarrow{\mathrm{\Omega }}\times {\overrightarrow{U}}_{ic}-\overrightarrow{\mathrm{\Omega }}\times \left(\overrightarrow{\mathrm{\Omega }}\times {\overrightarrow{r}}_{ic}\right)\right\},\hfill \end{array}$$

(9)

where $U_{ic}$, $U_r$, t, ${\rho }_a$, and ${\rho }_{ic}$ represent the velocity of the ice crystal, relative velocity between the surrounding flow and ice crystal, time, and densities of the air and ice crystals, respectively. An arrow above a variable indicates a vector component. ${\rho }_{ic}$ varies depending on the amount of water in the ice crystal during flight, defined as

$$\begin{array}{ccc}\hfill {\rho }_{ic}& =& \frac{W_i+W_w}{V_i+V_w},\hfill \end{array}$$

(10)

where $W_i$ represents the masses of ice, and $V_i$ and $V_w$ represent the volumes of ice and water, respectively. The second term on the right-hand side of Equation (9) represents the centrifugal and Coriolis forces, where Ω denotes the rotational speed of the coordinate system, and $r_{ic}$ denotes the radial position of the ice crystal. Further, the drag coefficient $C_D$ of the complete sphere [50] is denoted as

$$\begin{array}{ccc}\hfill C_D& =& \frac{24}{R{e}_{ic\_a}}\left(1+0.15R{e}_{ic\_a}^{0.687}\right),\hfill \end{array}$$

(11)

where $R{e}_{ic\_a}$ represents the particle Reynolds number based on $d_{ic}$, the viscosity of the surrounding air ${\mu }_a$, and $U_r$. When an ice crystal impinges on the walls and not the blade surface, it is reflected from the walls in a perfectly elastic collision. The time marching of the Equation (9) was computed using the Euler explicit method with the CFL number of $0.08$. The validity of this numerical procedure was confirmed in our previous studies [12,44,45].

4. Results and Discussion

4.1. Effects of the Ice Crystal Model, Relative Humidity, and MVD

The impinging ice crystal mass and water content are important for predicting ice crystal icing. This section evaluates the effects of these physical values on the blade per unit time, considering the effects of relative humidity and the initial MVD, as well as with and without the ice crystal model.

Table 3. Conditions of ice crystal trajectory and impingement simulation.

	Ice Crystal Model	Nusselt Number [-]	Relative Humidity [%]	MVD [μm]
Case 1-w/o	w/o	-	-
Case 2-$RH$40	w	$N{u}_F$	40
Case 2-$RH60$	w	$N{u}_F$	60	25, 50, 100, 200
Case 2-$RH80$	w	$N{u}_F$	80
Case 3-$RH60$	w	$N{u}_{WO}$	60

summarizes the classification of the present simulation cases. Case 1-w/o does not consider the state change of the ice crystal (i.e., no heat transfer to the ice crystal), whereas Cases 2 and 3 consider the ice crystal model as the heat transfer estimations by $N{u}_F$ and $N{u}_{WO}$, respectively. The cases are indexed according to their conditions; for example, Case 2-$RH60$-$MVD25$ refers to the case when an ice crystal with an initial MVD of 25 μm flies through a flow field with a relative humidity of 60% while taking into account the ice crystal model and forced convection.

Table 4. Percentages of tracked particles impinging on the rotor blade.

	MVD 25 μm [%]	MVD 50 μm [%]	MVD 100 μm [%]	MVD 200 μm [%]
Case 1-w/o	$15.5$	$21.9$	$25.4$	$27.0$
Case 2-$RH$40	$15.6$	$22.0$	$25.4$	$27.0$
Case 2-$RH60$	$15.7$	$22.0$	$25.4$	$27.0$
Case 2-$RH80$	$15.8$	$22.0$	$25.4$	$27.0$
Case 3-$RH60$	$15.7$	$22.0$	$25.4$	$27.0$

lists the percentages of ice crystals that were judged to have been impinged on the rotor blade among the tracked particles. The percentage of impinging ice crystals increased with the increasing MVD because larger particles are less likely to avoid impingement on the blade due to inertia. Note that a significant difference was not confirmed in the number of impingements between cases in which the ice crystal size (MVD) was the same, indicating that the impingement number was dominantly affected by the initial ice crystal size.

The impinging ice crystal mass and water content under the effect of relative humidity are shown in Figure 5 for different MVDs. Figure 5a shows the case for an MVD of 25 μm. Compared with Case 1-w/o (without the ice crystal model), the total mass of the impinging ice crystals (summation of water and ice masses) becomes slightly larger with increasing relative humidity. However, the larger particles of more than 100 μm show almost the same values against the relative humidity. This is due to the increase in the ratio of the surface area to the volume, resulting in the heat input becoming more diminutive than the unit volume. Moreover, as the relative humidity increases, the ice mass decreases, whereas the water content increases. This is because an increase in humidity corresponds to an increase in the heat of condensation $Q_3$ from the vapor to the water film, resulting in a thickened water film. Figure 5b–d show the results for MVDs of 50–200 μm. The trends in relative humidity are similar to those for an MVD of 25 μm, but the water content significantly decreases in the case of larger ice crystals. Accordingly, in the icing simulation of an ice crystal, the relative humidity plays an important role in addition to MVD, especially for smaller MVD.

Figure 6 shows the relationship between relative humidity and water content for Case 2-$RH40$, Case 2-$RH60$, and Case 2-$RH80$. The water content decreases according to a power low as MVD increases, which can be used to estimate the MVD to obtain a completely melted state (i.e., water droplet). Present results represented a power-law of −1.82, although the index might be changed by the blade’s shape. In the following, we discuss the impingement characteristics of ice crystals for a relative humidity of 60%, which is a typical value.

Figure 7 shows the total mass of the impinged ice crystal on the blade and the water content at a relative humidity of 60% for different MVDs. The influence of the state model using Equation (5) (i.e., Case 3-$RH60$) is displayed by comparing the case “w/o model” (i.e., Case 1-w/o), as shown in Figure 7a. In Case 3-$RH60$-$MVD25$ and Case 3-$RH60$-$MVD50$, the sums of the mass of the ice and water are slightly larger than those in Case 1-w/o-$MVD25$ and Case 1-w/o-$MVD50$, respectively, because the atmospheric moisture surrounding the ice crystals condenses and is incorporated into the water film of the ice crystals. As the MVD increased, the total mass (i.e., the sum of water and ice masses) increased, and the difference between Case 1-w/o and Case 3-$RH60$ became negligible. However, the water content decreased as the size of the ice crystals increased, and the water content for Case 3-$RH60$-$MVD25$ was significantly higher than that for other ice crystal sizes. This is because the heat per unit mass of the ice crystal received from the flow field increases with the increasing surface area for smaller ice particles, increasing the mass of the water film per unit time.

Figure 7b shows the results of Case 1-w/o and Case 2-$RH60$ (i.e., the state model with Equation (7)). The mass of the impinging ice crystal and water content show similar trends to those in Figure 7a, whereas the water content is high for each MVD. This difference is responsible for the heat transfer; the Nusselt number is fixed at 2 in Equation (5), but it varies due to the flow field, in accordance with Equation (7). When the Nusselt number varied according to the flow field conditions (Figure 7b), the water content increased by $8.1$%, $4.1$%, $1.5$%, and $0.54$% for MVDs of 25 to 200, respectively, compared to the case where the Nusselt number was fixed at 2 (Figure 7a). Accordingly, we found that the Nusselt number affects the water content, especially for small ice crystals with an MVD of 25 μm. It is well known that water content is an important and sensitive parameter in the prediction of ice crystal accretion [28]. Therefore, in the following, the impingement characteristics of the ice crystals on the rotor blade surface will be discussed in detail, focusing on Case 2-$RH60$, where the forced convection is considered in accordance with actual phenomena.

4.2. Impingement Characteristics

Figure 8 compares the collection efficiency of the ice impingement on the pressure surface of the blade for different MVDs of 25, 50, 100, and 200 μm in Case 1-w/o and Case 2-$RH60$. The collection efficiency is defined as the impingement mass per unit time and area. We do not display the suction surface as ice crystal impingement rarely occurs here. Note that Figure 8 represents the impingement mass of the only ice part of the ice crystals covered with water film, as shown in Figure 1. As shown in Figure 8a, Case 2-$RH60$-$MVD25$ apparently affects the collection efficiency; the collection efficiency is large at the leading edge, tip side, and trailing edge of the hub side in Case 1-w/o-$MVD25$, whereas it is large at the leading edge in Case 2-$RH60$-$MVD25$. This is due to the melting of ice crystals during flight and a decrease in the mass of ice impingement. For an MVD of 50 μm (Figure 8b), the amplitude of the collection efficiency was different, but the distribution was similar. For MVDs of 100 and 200 μm (Figure 8c,d), the distributions were almost identical; this is because the amount of melted ice is less than that in the case of an MVD of 25 μm.

Figure 9a–d show the water content of the ice crystals by applying the state change model impinged with the pressure surface for MVDs of 25, 50, 100, and 200 μm in Case 2-$RH60$. The minimum and maximum values for the color contour range were set individually within a 0–90% span of the rotor blades to visualize the water content distribution because the water contents near the tip (90–100% span) exceed 90% in all cases. Furthermore, the water content deviation between the leading and trailing edges at the mid-span increased as the ice crystal size decreased (see color contour range). This is because the flight time of the ice crystals before impinging on the rotor blade is longer for ice crystals colliding with the trailing edge; that is, the melting time of the ice crystals is longer. However, by comparing with the same chord position, the water contents on the leading edge of the tip side and the trailing edge of the hub side are higher. This is due to the long flight time of the ice crystal, as discussed later.

Accordingly, the impingement of ice crystals tends to occur on the leading edge. Moreover, if the ice crystal is large, it also occurs at the trailing edge of the hub side. Near the tip, the water content is particularly high in the case of a small size, which may reduce the cooling rate inside the blade even with the same number of impingements; thus, the possibility of ice formation is lower than that in those positions. In the case of an MVD of 50 μm, the water content at the time of blade impingement is approximately 20%, which means that ice accretion is most likely to occur under these conditions [28]. In addition, from the above discussion, the water content of the impinging ice crystal varies significantly for small ice crystals of less than 100 μm (see Figure 7), indicating the importance of trajectory simulation before impingement. The following sections focus on the ice crystal impingement characteristics for the leading edge and trailing edge, respectively.

4.3. Impingement Characteristics at the Leading Edge of the Rotor Blade

Figure 10a,b show the time histories of the static temperature of the flow field around the flying ice crystals injected in the 50–99% span with the state change model for MVDs of 25 and 200 μm in Case 2-$RH60$. Here, we confirmed that the injection locations of the ice crystal were slightly different in the circumferential direction whereas the impinging locations on the blade were the same. In each case, the static temperature at the time of ice crystal injection was approximately 303 K, and there was no significant difference depending on the span position. The static temperature increased sharply immediately before impingement due to the influence of air around the stagnation point near the leading edge. Therefore, the flight time until the impingement dominantly affects the ice crystal content. In both cases, the static temperature of the ice crystal injected from the 99% span significantly fluctuated between $8.0$ and $9.0$ ms as pointed out by the red frame.

Figure 11a,b show the static temperatures of the flow field around the rotor blade and the ice crystal trajectories at MVDs of 25 μm (blue line) and 200 μm (red line) at the 50% and 99% spans, respectively. The ice crystal injected from the 99% span flowed near the suction surface, as shown by the red frame in Figure 11b; therefore, the static temperature decreased compared with that when injected from the 50% span. For the flight time until impingement, Figure 10 shows that ice crystals injected from the 99% span took 2.13 and 2.05 times longer than those injected from the 50% span for MVDs of 25 μm and 200 μm, respectively. Therefore, the increase in water content at the tip side shown in Figure 9 can be attributed to the longer flight time.

We discuss the cause of the change in flight time for ice crystals injected from the 99% span. Figure 12 shows the axial velocity distribution of the flow field at the position where the ice crystals are injected. The axial velocity of the 95% to 99% span was slower than the axial velocity of the mid-span because of the boundary layer near the casing of the axial fan. The axial velocity of the 99% span was $0.474$ times that of the mid-span, and the velocity in the circumferential direction was dominant in the flow field. When ice crystals were injected into the flow field, the initial velocity was set to be the same as the velocity of the flow field at the injection position, which increased the flying time of ice crystals through the 99% span. The same tendency was observed for other ice crystal sizes, and the higher water content on the tip side shown in Figure 9 was caused by the difference in the flight times of the ice crystals.

4.4. Impingement Characteristics at the Trailing Edge of the Rotor Blade

Figure 13 shows the time history of the static temperature of the flow field around the flying ice crystals injected from the 1 to 80% span for MVDs of 25 and 200 μm in Case 2-$RH60$. In both cases, the initial static temperature was approximately 303 K, regardless of the span position; however, after approximately 4 ms, the temperature increased, with the highest value injected from the 80% span. As shown in Figure 9a,d, the water content of the impinging ice crystal on the trailing edge barely changed in the span location, except near the tip and hub. For the flight time of the ice crystals before impingement on the rotor blades, the case in which the ice crystals were injected from a 1% span required the longest flight time for MVDs of 25 and 200 μm. Therefore, the higher water content of ice crystals impacting the trailing edge of the hub side can be attributed to the increased flight time, which is similar to the case described in the previous subsection.

Figure 14a shows the axial velocity distribution of the flow field near the mid-chord of the rotor blade represented in Figure 14b, and Figure 14b shows the relative velocity vectors and ice crystal trajectories for the MVDs of 25 and 200 μm in the 1% span. Figure 14a indicates that the axial velocity of the 1% span is 0.489 times that of the mid-span due to the influence of the boundary layer formed around the rotor blade hub. This leads to the fact that the flying time of the ice crystals injected from the 1% span is slightly longer than that of ice crystals injected from other span positions. A comparison of the ice crystal trajectories is shown in Figure 14b for MVD 25 μm, where the trajectory change following the flow field can be observed near the rotor blade, while for MVD 200 μm, the trajectory is relatively straight near the blade. This is because an MVD of 25 μm gives a smaller Reynolds number than an MVD of 200 μm, as indicated by Equation (11), and the drag coefficient increases, which strongly affects the ice crystal trajectories. Ice crystals injected from the 1% span are affected by the boundary layer only when they pass near the hub side of the rotor blade. Therefore, the increases in flight time and water content are smaller than those of ice crystals impinging on the tip side.

5. Conclusions

In this study, we developed an ice crystal trajectory simulation in a three-dimensional flow field around the rotor blade of an axial fan. Computations were performed for four cases of MVD and three cases of relative humidity, and the impingement characteristics were clarified by considering the state changes of the ice crystals during flight. The following conclusions were obtained:

• The ice crystal model adopted in this study can reasonably predict the time required for a frozen ice crystal to melt completely by considering natural convection.
• Forced convection in the flow field is an important factor in the prediction of ice crystal icing, as the water content at the time of rotor blade impingement increases by up to $8.1$% for an MVD of 25 μm compared to the cases where the Nusselt number is fixed at 2.
• The water content of ice crystals at the time of rotor blade impingement was found to vary significantly depending on the flight time. Around the mid-span, the water content increases from the leading edge to the trailing edge.
• At the same chord position, the leading edge of the tip side and trailing edge of the hub side showed a particularly high water content. This is attributed to the effect of the boundary layer created in the casing and hub of the axial fan, which increased the flight time compared to ice crystals flying in other positions.
• More condensation heat was generated with an increase in the relative humidity, which increased the amount of water vapor taken into the ice crystal and the amount of ice melted, which, in turn, increased the water content of the ice crystal at the time of impingement.

The tendencies of the latter three conclusions became more pronounced as the ice crystal size decreased. If the ice crystal size is particularly large, ice crystal icing is most likely to occur at the leading edge of the rotor blade and the trailing edge of the hub side. In the case of MVD 50 μm, it was easier to adhere to the rotor blade, except for near the tip region because the water content at the tip side shows particularly high values. The results of this study demonstrate that the impingement characteristics of ice crystals vary greatly depending on the span and chord position, and it is essential to consider all three dimensions in icing simulations.