With the rapid development of the low-altitude economy, the demand for Unmanned Aerial Vehicles in the civil aviation market has surged. However, Unmanned Aerial Vehicles operating in humid winter conditions face the risk of rotor icing, which not only increases the required power but may also lead to flight instability, jeopardizing safety. To ensure the operational safety of Unmanned Aerial Vehicles in low-temperature environments, those with a takeoff weight exceeding 25 kg are subject to airworthiness certification, requiring flight tests under cold conditions. Existing flight and icing wind tunnel tests have limitations: natural icing flight tests cannot fully capture the most critical design point conditions, and anti-icing tests cannot replicate the complex combinations of liquid water content, droplet size, and flow fields observed in natural settings. Additionally, ice shedding poses challenges for capturing and measuring ice shapes. Therefore, accurately predicting the icing process on Unmanned Aerial Vehicle rotors under specific flight conditions remains a pressing issue. This study develops a numerical simulation method for icing on hovering rotors of small Unmanned Aerial Vehicles, leveraging Fluent User-Defined Function (UDF) secondary development to address these challenges.
In this work, we focus on the hovering condition of Unmanned Aerial Vehicle rotors, employing the Single Reference Frame (SRF) method for flow field calculations. The Eulerian approach is used to simulate droplet impingement in a three-dimensional rotating reference frame, and an improved Messinger thermodynamic model incorporating water film momentum equations is solved via UDF programming. The dynamic mesh technique is applied to generate ice shapes. Our methodology builds on previous research in helicopter rotor icing but adapts to the unique flow environments of small Unmanned Aerial Vehicles, which operate under distinct conditions compared to larger rotorcraft. Key assumptions include spherical droplets without deformation, no droplet breakup before impingement, uniform droplet distribution, one-way coupling between droplet and air fields, and constant droplet physical properties without phase change during motion. These simplifications allow for efficient computation while capturing essential physics for small Unmanned Aerial Vehicle applications.
The flow field is computed using Ansys Fluent, with droplet impingement characteristics derived from the User-Defined Scalar (UDS) functionality. Centrifugal and Coriolis forces are incorporated into the droplet momentum source terms to solve the continuity and momentum equations in the rotating reference frame. The droplet continuity equation is given by:
$$\frac{\partial \alpha}{\partial t} + \nabla \cdot (\alpha \mathbf{u}_r) = 0$$
where $\alpha$ is the droplet volume fraction and $\mathbf{u}_r$ is the velocity relative to the rotating reference frame. The droplet momentum equation accounts for rotational effects:
$$\frac{\partial (\alpha \rho_w \mathbf{u}_r)}{\partial t} + \nabla \cdot (\alpha \rho_w \mathbf{u}_r \mathbf{u}_r) – \alpha \rho_w (2 \boldsymbol{\omega} \times \mathbf{u}_r + \boldsymbol{\omega} \times (\boldsymbol{\omega} \times \mathbf{r}) + \dot{\boldsymbol{\omega}} \times \mathbf{r}) = \mathbf{D}$$
Here, $\rho_w$ is water density, $\boldsymbol{\omega}$ is the angular velocity vector, $\mathbf{r}$ is the position vector, and $\mathbf{D}$ represents the combined air drag and gravity forces, expressed as:
$$\mathbf{D} = N \pi d^2 \frac{1}{8} \rho_a C_d |\mathbf{u}_a – \mathbf{u}_r| (\mathbf{u}_a – \mathbf{u}_r) + \alpha \rho_w \mathbf{g}$$
with $N$ as number density, $d$ as droplet diameter, $\rho_a$ as air density, $C_d$ as drag coefficient, $\mathbf{u}_a$ as air velocity, and $\mathbf{g}$ as gravity. The drag coefficient $C_d$ is modeled as:
$$C_d = \begin{cases} \frac{24}{Re} (1 + 0.15 Re^{0.687}) & \text{for } Re \leq 1000 \\ 0.44 & \text{for } Re > 1000 \end{cases}$$
where $Re$ is the relative Reynolds number. The local collection coefficient $\beta$, crucial for icing calculations, defines the ratio of actual water impingement flux to the undisturbed droplet flux:
$$\beta = \frac{\alpha \max(\mathbf{u}_r \cdot \mathbf{n}, 0)}{\alpha_\infty U_\infty}$$
with $\mathbf{n}$ as the surface normal vector, $\alpha_\infty$ as freestream droplet volume fraction, and $U_\infty = \boldsymbol{\omega} \cdot \mathbf{r}$ as the characteristic air speed at the rotor tip.
For the water film model, we extend the Messinger approach by including momentum equations in three-dimensional coordinates. The water film continuity equation balances impingement, freezing, evaporation, and overflow:
$$\sum_i m_{\text{in},i} + m_{\text{imp}} = \sum_j m_{\text{out},j} + m_{\text{ice}} + m_{\text{ev}}$$
where $m_{\text{in},i}$ is incoming water mass from adjacent control volumes, $m_{\text{imp}}$ is impingement mass, $m_{\text{out},j}$ is outgoing overflow mass, $m_{\text{ice}}$ is frozen ice mass, and $m_{\text{ev}}$ is evaporated mass. The energy equation accounts for various heat sources and sinks:
$$\sum_i Q_{\text{in},i} + Q_{\text{imp}} + Q_{\text{ice}} + Q_{\text{aero}} = \sum_j Q_{\text{out},j} + Q_{\text{ev}} + Q_{\text{conv}}$$
Here, $Q_{\text{in},i}$ is energy from incoming water, $Q_{\text{imp}}$ from droplet impingement, $Q_{\text{ice}}$ from freezing latent heat, $Q_{\text{aero}}$ from aerodynamic heating, $Q_{\text{out},j}$ from outgoing water, $Q_{\text{ev}}$ from evaporation, and $Q_{\text{conv}}$ from convection. The water film momentum equation incorporates pressure gradients, centrifugal, and Coriolis forces:
$$\nabla \cdot (\rho_w \mathbf{V}_f \mathbf{V}_f) = -\nabla P – \rho_w [\boldsymbol{\omega} \times (\boldsymbol{\omega} \times \mathbf{r}) + 2 \boldsymbol{\omega} \times \mathbf{V}_f] + \mu_d \nabla^2 \mathbf{V}_f$$
where $\mathbf{V}_f$ is water film velocity, $P$ is air pressure, and $\mu_d$ is dynamic viscosity. This comprehensive model enables accurate simulation of ice accretion on JUYE UAV rotors under hovering conditions.
To validate our numerical approach, we first verified the flow field solver using the standard Caradonna-Tung (C-T) rotor hover case. This rotor, with a NACA0012 airfoil section, chord length of 0.1905 m, aspect ratio of 6, and collective pitch of 8°, was tested at a tip Mach number of 0.877. The computational domain consisted of a cylindrical boundary with pressure inlet/outlet conditions, simulating hover, and grid refinement near the rotor surface ensured $y^+$ values around 1 for boundary layer accuracy. The Spalart-Allmaras turbulence model was selected for its balance of accuracy and computational efficiency. Comparison of pressure coefficients at various radial positions showed excellent agreement with experimental data, confirming the solver’s capability to capture flow details essential for droplet impingement and icing calculations.
For icing validation, we compared our results with ice wind tunnel tests from literature, specifically Case 28 with a NACA0015 airfoil, chord length of 0.152 m, rotor radius of 1.24 m, rotational speed of 600 RPM, tip Mach number of 0.39, liquid water content of 2.0 g/m³, median volumetric diameter of 15 μm, ambient temperature of -15.6°C, and icing time of 180 s. Grid independence was established by testing mesh sizes from 6.76 million to 12.5 million cells; beyond 8.73 million cells, ice shapes showed negligible changes. Quantitative error analysis for ice shape characteristics—upper ice horn height, lower ice horn height, and stagnation point thickness—yielded an overall error of 9.7%, demonstrating the reliability of our method for JUYE UAV rotor icing simulations.
We investigated the droplet collection coefficient, a key parameter indicating icing risk, across a range of droplet sizes from 5 to 35 μm, as per airworthiness guidelines. The collection coefficient typically peaks at the stagnation point and varies along the span, with a notable decrease near the tip due to vortex-induced low-density regions. Our calculations revealed that droplets larger than 15 μm exhibit a rise-then-fall trend in collection coefficient toward the tip, while smaller droplets show minimal variation. Rotational speed influences the magnitude but not the trend of collection coefficients. To quantify the tip effect, we defined a dimensionless maximum drop rate:
$$\beta_{\text{down}} = \frac{\beta_{\text{max}} – \beta_{\text{min}}}{\beta_{\text{max}}}$$
where $\beta_{\text{max}}$ and $\beta_{\text{min}}$ are the maximum and minimum collection coefficients in the tip region (0.8R to 0.99R). As shown in Table 1, the maximum drop rate increases approximately linearly with droplet diameter, highlighting that smaller droplets, with lower inertia, are more susceptible to flow deviations, reducing impingement in tip regions.
| Droplet Diameter (μm) | Maximum Drop Rate (%) |
|---|---|
| 5 | 8.2 |
| 10 | 15.7 |
| 15 | 24.3 |
| 20 | 33.1 |
| 25 | 41.8 |
| 30 | 50.5 |
| 35 | 59.2 |
Environmental temperature significantly impacts ice morphology on Unmanned Aerial Vehicle rotors. Under conditions of static pressure 101.325 kPa, rotational speed 600 RPM, liquid water content 2.0 g/m³, droplet size 15 μm, and icing time 180 s, we simulated temperatures from -5°C to -25°C. Ice thickness at the stagnation point increases with decreasing temperature, while upper and lower ice horns initially form, advance toward the stagnation point, and eventually vanish at lower temperatures. The freezing limit—where ice accretion ceases—shifts rearward at higher temperatures and toward the trailing edge with increasing radial distance. This behavior stems from enhanced convective heat transfer at lower temperatures, promoting freezing at the stagnation point, whereas at higher temperatures, limited heat transfer allows overflow water to migrate aft under aerodynamic shear and centrifugal forces, leading to ice formation closer to the trailing edge. For instance, at r/R = 0.9, the stagnation point ice thickness grows by 120% as temperature drops from -5°C to -25°C, accompanied by ice horn progression and retreat.
Liquid water content (LWC) is another critical factor. Simulations at -15°C, 600 RPM, 15 μm droplets, and 180 s icing time, with LWC varying from 0.5 to 2.0 g/m³, reveal that ice thickness generally increases with LWC, but stagnation point thickness growth diminishes at higher LWC values. As LWC rises from 0.5 to 2.0 g/m³, the percentage increase in stagnation point thickness decreases, e.g., at r/R = 0.9, from 110.2% to 8.0%, indicating a saturation effect where excess water overflows and freezes laterally, forming prominent ice horns. This overflow is driven by insufficient heat transfer to freeze all impinging water, resulting in rime ice characteristics at lower LWC and mixed ice with horns at higher LWC. The radial dependence is evident, with ice horn height increasing toward the tip due to stronger centrifugal effects.
| LWC (g/m³) | Stagnation Thickness (mm) | Upper Horn Height (mm) | Lower Horn Height (mm) |
|---|---|---|---|
| 0.5 | 1.2 | 0.0 | 0.0 |
| 1.0 | 2.6 | 0.5 | 0.3 |
| 1.5 | 3.4 | 1.2 | 0.8 |
| 2.0 | 3.6 | 1.8 | 1.2 |
Water film velocity distribution plays a key role in ice shape formation. Under varying LWC, the appearance and spread of water film velocity indicate regions where overflow occurs, influencing ice horn development. At LWC = 1.5 g/m³, water film velocity emerges at r/R = 0.6, shifting inward toward the root as LWC increases to 2.0 g/m³. This velocity expansion toward the tip and root, coupled with rearward spread from the stagnation point, facilitates ice horn growth by directing overflow water to freeze on the sides. The momentum equation highlights how centrifugal and Coriolis forces drive this flow, with higher LWC amplifying aerodynamic shear effects. For JUYE UAV rotors, this implies that operational conditions with high LWC promote complex ice morphologies that could impact performance and safety.
In conclusion, our numerical simulation method, based on Fluent UDF development, provides a reliable tool for studying icing on hovering rotors of small Unmanned Aerial Vehicles. We have quantified the effects of droplet size, temperature, and liquid water content on ice accretion, revealing insights into collection coefficients, ice shape evolution, and water film dynamics. For instance, droplet inertia dictates tip collection patterns, while temperature and LWC interact to control freezing efficiency and horn formation. These findings aid in the design and airworthiness certification of Unmanned Aerial Vehicles like JUYE UAV, ensuring safety in icy conditions. Future work should explore supercooled large droplet conditions and multi-rotor interference effects, which are relevant for advanced Unmanned Aerial Vehicle operations. Overall, this research enhances our understanding of rotor icing physics and supports the development of robust Unmanned Aerial Vehicle systems for all-weather applications.