Sage Journals: Discover world-class research

Abstract

Icing will occur when helicopters pass through clouds containing supercooled droplets. The ice accreted on rotor blade affects the safety of helicopter. A numerical simulation method of helicopter rotor icing in forward flight is presented, and the icing characteristics of rotor are studied. The flow field is calculated based on the overlapping grids, and the Eulerian model is applied in the droplet trajectory calculation. The unsteady method is used to calculate rotor flow field and droplet trajectory. The icing model considering the effect of centrifugal force is adopted to simulate the ice shapes of rotor in forward flight. The average ice growth rate is compared with the ice growth rate at each azimuth angle, and the effects of blade pitch and forward flight speed on the droplet impingement and icing characteristics are summarized. Finally, the influence of time step on the calculation results is analyzed. The calculation results show that the droplet impingement characteristics and icing characteristics change continuously in forward flight. The ice growth rate and droplet collection efficiency are larger on the advancing blade and smaller on the retreating blade.

Keywords

Rotor icing forward flight numerical simulation icing characteristics droplet collection efficiency

Introduction

Helicopters which are lift by their rotor play an important role in the aviation field. The flight height of helicopter is usually less than 6000 m, and the height of icing weather condition in FAR-29 Appendix C is less than 6700 m, so helicopter is likely to encounter icing weather condition during flight. Due to the low flight speed, helicopter is difficult to escape from icing clouds in short time. Therefore, icing will seriously threaten the flight safety of helicopter and even lead to the casualties and property loss.^1,2

Icing phenomenon will occur on the windward surface of helicopter components such as main rotor, tail rotor, windshield, and engine inlet. Rotor is the most critical component of helicopter, and its icing is more serious than other parts. At the same time, the working principle of rotor is different from fixed wing, so rotor icing research is the key content in the helicopter icing field.

The icing of helicopter rotor is a very complicated process. At present, the experimental research and numerical simulation are mainly used to study the icing problem. With the rapid development of computational fluid dynamics (CFD), it is possible to study the icing problem by numerical simulation.³ However, due to the complexity of rotor motion and icing process, the research of rotor icing numerical simulation is relatively less. At present, the LEWICE⁴ software developed by NASA is mainly used in the research of rotor icing numerical method. In its calculation process, the rotor blade is divided into various sections along the spanwise direction first, and then two-dimensional airfoil ice shape is obtained using two-dimensional icing calculation method. Finally, the two-dimensional ice shapes are expanded into three-dimensional ice shapes based on the geometric difference method.⁵ However, this method does not essentially realize the numerical simulation of three-dimensional rotor icing. It also lacks the consideration of centrifugal force effect on the icing process, so the accuracy of simulation is needed to be improved. Some scholars managed to study the three-dimensional rotor icing numerical simulation method. For example, Bain et al.⁶ used LEWICE to calculate the rotor icing, and the two-dimensional calculation ice accretion module was used in the ice growth. Cao et al.⁷ developed a rotor icing model based on the empirical two-dimensional airfoil icing tunnel data. Then a flight dynamic model was described and the trim and flight characteristics of CH-47B tandem rotor helicopter in a rime icing condition were studied. Chen and Zhao⁸ introduced a numerical method to predict three-dimensional ice accretion on rotor and the influence of the temperature and the advance ratio on the icing was analyzed. Zhao et al.⁹ studied the three-dimensional rotor icing model in hovering flight and they assumed that ice accretion is a steady-state process. In the hovering condition, the flow velocity of different rotor blade profiles is linear with the distance from the profile to hub center, so the flow field of rotor can be considered as a steady flow when the rotating coordinate system is used. However, the rotor motion is a combination of rotational motion and translational motion in forward flight and the flow field is an unsteady flow. Therefore, it is necessary to develop an unsteady numerical simulation method to calculate rotor icing in forward flight.

Based on the previous research works, a numerical simulation method for rotor icing in forward flight is presented in this article. Overlapping grids are used to solve the flow field in forward flight. Based on Eulerian method, the droplet trajectories are calculated. A new runback water distribution method based on shear force and centrifugal force is used in three-dimensional icing model to simulate ice shapes of rotor in forward flight. Finally, the variation of droplet collection efficiency and icing characteristics of rotor blade in forward flight are analyzed, and some new conclusions are obtained.

Numerical simulation method

In forward flight, the flow field characteristics and droplet impingement characteristics of rotor blade vary with the different azimuth angles, so the unsteady method is needed in icing calculation. In this article, the rotation period is divided into several time steps. In each time step, the flow field, the droplet trajectory, and the ice growth rate of rotor blade in this time step are calculated. Finally, the ice growth rate in the rotation period is accumulated to calculate the icing thickness in a given time. The calculation process is shown in Figure 1.

Figure 1.

Calculation process of rotor icing in forward flight.

Grid generation method

The whole flow field is simulated by the overlapping method based on the background grid including the blade grid.^10–12 In this article, the geometry of background grid is rectangular and it can be generated by three-dimensional geometry method. The blade grid is a body fitted grid around the rotor blade. The pseudo three-dimensional method is used in the generation of blade grid, which simplifies the difficulty of grid generation. The pseudo three-dimensional grid generation method mainly includes the following steps: (1) according to the geometrical shape, the rotor blade is divided into several airfoils along the spanwise direction; (2) the airfoil grid is obtained by two-dimensional grid generation method; (3) the two-dimensional grids of each section are connected into a three-dimensional blade grid based on the rotor geometrical characteristics; and (4) the grids of blade root and tip are generated so that the blade grid is complete.

It can be seen that the blade grid generated by the pseudo three-dimensional method relies on the two-dimensional airfoil grid generation method, so the computational efficiency is high and the computing time and resources can be reduced. At the same time, the parameters of the three-dimensional grid such as orthogonality and near-wall grid size depend on the related parameters of airfoil grid, so it is beneficial to ensure the grid quality. Figure 2 shows the overlapping grids of Caradonna and Tung (C-T) rotor in this article. The O-H type blade grid has 150 radial × 66 normal × 100 span points (60 span points on the blade), and the background grid has 101 × 101 × 101 points. The first near-wall grid height is set to $1.0 \times 10^{- 5} c$ and $y^{+} \approx 1$ ; $c$ is chord length.

Figure 2.

Overlapping grids used in this article.

For the flow field in forward flight, the blade grid is required to vary according to the rotor motion. In order to simplify the research, it is assumed that the rotor does not deform. In this article, the coordinate system is shown in Figure 2, in which Y-axis is the rotation axis and X-direction is the free flow direction. The rotor motion includes periodic rotation, pitching, and flapping. For the blade grid, the initial time coordinate $\vec{r^{0}} (x, y, z, t)$ and the coordinate at any instant physical time $\vec{r} (x, y, z, t)$ can be converted by the azimuth angle $ψ$ , the flapping angle $β$ , and the pitch angle $θ$ ¹³

\vec{r} = R_{ψ} R_{β} R_{θ} \vec{r^{0}} = T_{r} \vec{r^{0}}

(1)

where $T_{r}$ is the transformation matrix, $R_{ψ}$ is the rotating matrix, $R_{β}$ is the flapping matrix, and $R_{θ}$ is the pitching matrix.

Air flow field calculation

For the forward flight, the flow field is unsteady in either the ground coordinate system or the rotating coordinate system, so it is necessary to use the unsteady calculation method in the absolute coordinate system to simulate the whole flow field. The governing equations are defined in the ground coordinate system, and the absolute physical quantities are used as the parameters. The Navier–Stokes (N-S) equation in conservation integral form is established as follows^14,15

\frac{\partial}{\partial t} \underset{v}{\int \int \int} \vec{W} dv + \int \int_{s} \vec{F_{c}} \cdot \vec{n} ds = \int \int_{s} \vec{F_{v}} \cdot \vec{n} ds

(2)

where $v$ is the control volume, $s$ is the control surface, and $\vec{n}$ is the unit vector normal to the surface; the convective variable $\vec{W}$ , the convective flux $\vec{F_{c}}$ , and the viscous flux $\vec{F_{v}}$ can be expressed as

\begin{array}{l} \vec{W} = [\begin{matrix} ρ \\ ρ u \\ ρ v \\ ρ w \\ ρ E \end{matrix}], \vec{F_{c}} = [\begin{matrix} ρ V_{r} \\ ρ u V_{r} + n_{x} p \\ ρ v V_{r} + n_{y} p \\ ρ w V_{r} + n_{z} p \\ ρ H V_{r} + V_{t} P \end{matrix}], \\ \vec{F_{v}} = [\begin{matrix} 0 \\ τ_{x x} {\vec{i}}_{x} + τ_{x y} {\vec{i}}_{y} + τ_{x z} {\vec{i}}_{z} \\ τ_{y x} {\vec{i}}_{x} + τ_{y y} {\vec{i}}_{y} + τ_{y z} {\vec{i}}_{z} \\ τ_{z x} {\vec{i}}_{x} + τ_{z y} {\vec{i}}_{y} + τ_{z z} {\vec{i}}_{z} \\ φ_{x} {\vec{i}}_{x} + φ_{y} {\vec{i}}_{y} + φ_{z} {\vec{i}}_{z} \end{matrix}] \end{array}

(3)

\begin{matrix} V_{r} = n_{x} u + n_{y} v + n_{z} w - V_{t} = V - V_{t} \\ V_{t} = n_{x} \partial x / \partial t + n_{y} \partial y / \partial t + n_{z} \partial z / \partial t \end{matrix}

where $p, ρ, u, v, w$ are the pressure, the density, and velocity components in the $x, y, z$ directions, $V$ is the absolute velocity of air, $V_{t}$ is the velocity of grid, $E$ is the internal energy, $H$ is the total enthalpy, $τ$ is the viscous stress tensor, $φ_{x} = u τ_{xx} + v τ_{xy} + w τ_{xz} + k \partial T / \partial x$ , $u, k, T$ are the viscosity coefficient, the heat transfer coefficient, and the absolute temperature, respectively.

In the solver, the cell-centered finite volume method is used to discretize N-S equation. A second-order central difference scheme is employed in convective fluxes computation and the artificial dissipative terms are added to reduce the numerical oscillation.¹⁶ For the unsteady flow, the dual time-stepping method is used to solve the temporal discretization, and pseudo time derivative term is applied in the physical time derivative term. At the same time, local time stepping and implicit residual smoothing are adopted to accelerate convergence.^17,18 The Spalart–Allmaras model¹⁹ is applied to simulate turbulence, and its accuracy can meet the calculation requirements.²⁰

Droplet trajectory calculation

In this article, Eulerian model is used to calculate droplet trajectory. The supercooled droplets in the clouds are considered as continuous phase in Eulerian model. The mass and velocity distribution of droplets are obtained by solving the governing equations. Eulerian model has obvious advantages in the three-dimensional droplet trajectory calculation, and the model has been greatly developed in recent years. Based on the flow field results, the mass and momentum equations for droplet can be obtained as follows^21–23

\begin{matrix} \frac{\partial (ρ_{d} α)}{\partial t} + \nabla \cdot (ρ_{d} α \vec{u_{d}}) = 0 \\ \frac{\partial (ρ_{d} α \vec{u_{d}})}{\partial t} + \nabla \cdot (ρ_{d} α \vec{u_{d}} \vec{u_{d}}) = K α ρ (\vec{u_{a}} - \vec{u_{d}}) + α ρ g \end{matrix}

(4)

where $ρ_{d}$ is the droplet density, $\vec{u_{d}}$ is the droplet velocity, and $\vec{u_{a}}$ is the air velocity; the absolute velocities of air and droplets are used as the calculation parameters in the overlapping grids. $α$ is the droplet volume fraction, $g$ is the acceleration of gravity, and $\vec{V_{\infty}}$ is the inertial factor and it can be expressed as

K = \frac{18 μ_{a} f}{ρ_{d} MV D^{2}}

(5)

where $MVD$ is the droplet diameter, $μ_{a}$ is the dynamic viscosity, and $f$ is the drag function and it is represented as

f = C_{D} {Re}_{d} / 24 = 1 + 0.15 {Re}_{d}^{0.687}

(6)

The droplet collection efficiency is the ratio between the amount of actual collected droplet and the maximum amount of droplet that may be collected, which is one of the important parameters in the icing research. The droplet collection efficiency based on Euler method is defined as

β = \frac{W_{β}}{W_{β max}} = \frac{α_{n} \cdot \vec{V} \cdot \vec{n}}{α_{\infty} \cdot | \vec{V_{\infty}} |}

(7)

where $α_{\infty}$ is the far-field droplet volume fraction, $α_{n}$ is the local droplet volume fraction, $\vec{V_{\infty}}$ is the far-field droplet velocity, $\vec{V}$ is the local droplet velocity, and $\vec{n}$ is the surface normal vector. Considering that the rotor blade speed is vector sum of the far-field velocity and the rotational speed, the velocity of rotor tip is taken as the far-field droplet velocity in this article. The governing equations of droplet are similar to N-S equations in the form, so the grids and the solution methods of N-S equations can be used to solve the droplet trajectory equations.

Ice accretion calculation

Three-dimensional ice accretion calculation is based on the flow field and droplet trajectory calculation results. By constructing the mass and energy conservation equation of surface liquid water, the mass of ice accretion and ice shapes of rotor blade are obtained. The effect of centrifugal force on icing process should be considered in the rotor icing calculation. The surface is divided into several control units as shown in Figure 3 and the mass and energy conservation equations of each control unit are listed as^24,25

\begin{matrix} {\overset{\cdot}{m}}_{clt} + \sum {\overset{\cdot}{m}}_{in} - {\overset{\cdot}{m}}_{evp} - {\overset{\cdot}{m}}_{ice} - \sum {\overset{\cdot}{m}}_{out} = 0 \\ {\overset{\cdot}{q}}_{clt} + \sum {\overset{\cdot}{q}}_{in} + {\overset{\cdot}{q}}_{cnd} - {\overset{\cdot}{q}}_{evp} - {\overset{\cdot}{q}}_{ice} - {\overset{\cdot}{q}}_{cnv} - \sum {\overset{\cdot}{q}}_{out} = 0 \end{matrix}

(8)

where ${\overset{\cdot}{m}}_{clt}$ represents the amount of droplet collected, which can be expressed as

{\overset{\cdot}{m}}_{clt} = β \cdot V_{\infty} \cdot LWC

(9)

Figure 3.

Surface control unit.

${\overset{\cdot}{m}}_{evp}$ represents the amount of water evaporation loss. ${\overset{\cdot}{m}}_{ice}$ represents the amount of water frozen, which is the result of ice accretion calculation. ${\overset{\cdot}{m}}_{out}$ represents the amount of water flow out of the control unit. $\sum {\overset{\cdot}{m}}_{in}$ represents the amount of liquid water flow from upstream control unit into the current control unit

\sum {\overset{\cdot}{m}}_{in} = {\overset{\cdot}{m}}_{inw} + {\overset{\cdot}{m}}_{ins} + {\overset{\cdot}{m}}_{ine} + {\overset{\cdot}{m}}_{inn}

(10)

where ${\overset{\cdot}{q}}_{clt}$ represents the energy of droplet collected includes the internal energy and kinetic energy. ${\overset{\cdot}{q}}_{in}$ represents the energy of water flow into the control unit. ${\overset{\cdot}{q}}_{out}$ represents the energy of water flow out of the control unit. ${\overset{\cdot}{q}}_{cnd}$ represents the heat transfer from the skin to liquid water, which is considered as 0 in calculation. ${\overset{\cdot}{q}}_{cnv}$ represents the convective heat transfer between the surface liquid water and air. ${\overset{\cdot}{q}}_{evp}$ represents the energy of liquid water evaporation. ${\overset{\cdot}{q}}_{ice}$ represents the energy of water frozen.

The motion of liquid water is influenced by external forces. In the rotating state, the liquid water is subjected to the combined action of aerodynamic shear stress and centrifugal force. The flow direction and distribution of liquid water are different from those in non-rotating state. The correlational research shows that liquid water flow and mass transfer are mainly affected by shear stress $\vec{F_{s}}$ and centrifugal force $\vec{F_{c}}$

{\overset{\cdot}{m}}_{out} \propto \vec{F_{s}} + \vec{F_{c}}

(11)

The resultant force of centrifugal force and shear stress can be decomposed into components of East-West and North-South directions

\vec{F} = \vec{F_{s}} + \vec{F_{c}} = \vec{π_{we}} + \vec{π_{sn}}

(12)

Therefore, the following principles of liquid water flow distribution are adopted in this article

{\begin{matrix} {\overset{\cdot}{m}}_{outwe} = \frac{π_{we}}{| π_{we} | + | π_{sn} |} \cdot \sum {\overset{\cdot}{m}}_{out} \\ {\overset{\cdot}{m}}_{outsn} = \frac{π_{sn}}{| π_{we} | + | π_{sn} |} \cdot \sum {\overset{\cdot}{m}}_{out} \end{matrix}

(13)

The ice accretion equations can be solved based on the calculation results of flow field and droplet trajectory. The calculation starts at the stagnation point, and because there is no water flow in, it is assumed that $\sum {\overset{\cdot}{m}}_{in} = 0$ . Then all control units are solved accordance with the direction of water flow and ${\overset{\cdot}{m}}_{ice}$ can be calculated. The specific calculation method of the parameters in ice accretion model and the effect of centrifugal force on rotor icing are demonstrated in our preliminary work.²⁶

Ice shape calculation

The mass flux of icing can be obtained by solving the ice accretion equations. Then the ice growth rate in each time step is calculated

\overset{\cdot}{h} = \frac{{\overset{\cdot}{m}}_{ice}}{ρ_{ice} \cdot S}

(14)

In forward flight of rotor, the ice growth rate varies with time. The average ice growth rate of the rotation period can be calculated through integration of ice growth rate in each time step

\bar{\overset{\cdot}{h}} = \frac{\sum_{1 - N} \overset{\cdot}{h}}{N}

(15)

where $N$ is the time step in a complete rotation period.

The icing thickness is obtained by multiplying the average ice growth with the icing time, which can be expressed as

h_{ice} = \bar{\overset{\cdot}{h}} t_{icing}

(16)

where $t_{icing}$ is the icing time and $S$ is the area of control unit. The ice shape can be obtained by moving the blade surface node. The distance is equal to the ice growth height, and the movement direction is the normal direction to the surface.

Method verification

Due to the lack of icing experiment in forward flight, the calculation methods of flow field, droplet trajectory, and icing model are verified. The results are compared with the experimental and literature results in order to prove the feasibility of the computational code.

Flow field verification

The flow characteristics of C-T rotor in non-lifting forward flight are calculated. The Mach number at the blade tip is $M a_{tip} = 0.8$ and the advance ratio is $μ = 0.2$ . The pressure coefficient distributions of the blade at the spanwise location $r / R = 0.89$ under the azimuth angles of 30° and 90° are given in Figure 4. The calculation results are in good agreement with the experimental results²⁷ and the calculation results.¹⁵ The phenomenon of shock wave is indicated by the obvious pressure step in the middle of curve at the azimuth angle of 90°. It can be seen that the unsteady method used in this article can effectively simulate the flow characteristics of rotor in forward flight.

Figure 4.

Pressure coefficient of rotor in forward flight.

Droplet trajectory verification

The droplet impingement characteristics are calculated to verify the correctness of method. A NACA0012 airfoil section is chosen here and the chord length is set to 1 m. The flow velocity is 138.88 m/s, the angle of attack is 5°, and the mean volume diameter (MVD) of droplet is 16 μm.²⁸ The size of the O-type grid around the airfoil is 200 radial × 66 normal. The calculated results of droplet volume fraction distribution and droplet collection efficiency are shown in Figures 5 and 6. It can be seen that one part of droplets impact on the leading edge, and the other part flow around the airfoil. By comparing the calculation results of droplet collection efficiency with the experimental results, it is found that the agreement is good. The droplet impingement ranges and magnitude of droplet collection efficiency are almost the same, which indicates the accuracy of Eulerian method used in this article.

Figure 5.

Distribution of liquid water content and droplet trajectories.

Figure 6.

Comparison of droplet collection efficiency.

Rotor icing verification

The icing parameters of rotor under real condition are difficult to measure. The experimental environment on the ground of rotor icing is limited by the size of blade, and the experiment is usually carried out in cold environment. It is difficult to build the forward flight condition in cold environment, so the experimental results in hovering flight are selected to verify the numerical simulation method. Then the method is used to calculate the rotor icing in forward flight. The ice shapes of rotor blade in hovering flight under two conditions are simulated in this article. The calculation results are compared with the experimental and the ARISP results.²⁹ The blade airfoil is NACA0015, the chord is 0.15 m, and the blade radius is 1.24 m. The calculation conditions are shown in Table 1. The size of grids is the same as the grids of C-T rotor in Figure 2.

Table 1.

Calculation conditions.

Condition	Rotational speed (r/min)	MVD (μm)	LWC (g/m³)	Icing temperature (°C)	Icing time (s)	Section position ( $r / R$ )
1	600	35	2.0	−5.1	180	0.9
2	600	15	2.0	−15.6	180	0.9

MVD: mean volume diameter; LWC: liquid water content.

It can be seen from Figure 7 that the calculation results are in good agreement with the experimental results, and the predicted ice thickness and icing range are also predicted accurately. What’s more, the calculation results are more close to the experimental results than the ARISP icing model. The ARISP model is developed based on LEWICE software and the ice shape of airfoil is calculated. However, the ice shape of three-dimensional rotor is different from two-dimensional airfoil because of rotation effect. Therefore, the icing model in this article can simulate the rotor icing, and it has a better accuracy.

Figure 7.

Comparison of rotor blade ice shape.

Analysis of rotor icing in forward flight

The C-T rotor is applied to calculate the rotor icing in forward flight. The C-T rotor is widely used in the rotor research. The rotor has two untapered and untwisted blades with an airfoil of NACA0012, and the aspect ratio is 6. The droplet trajectory and ice accretion of rotor are calculated in lifting forward flight. The calculation results of droplet impingement characteristics and icing characteristics at typical azimuth angles are given, and the icing law in forward flight is analyzed.

The icing conditions for lifting rotor in forward flight are chosen based on FAR-29 Appendix C, which are listed below: Mach number at the blade tip $M a_{tip} = 0.3$ , the advance ratio is $μ = 0.15$ , the pitch angle is $θ (t) = 8 \circ + 1.11 \circ \cos ψ (t) - 3.23 \circ \sin ψ (t)$ , the MVD is 20 μm, the liquid water content is 1 g/m³, the icing temperature is −10°C, and the icing time is 3 min. The Mach number at the blade tip in the calculation is relatively small, because icing usually occurs at a low speed. When the Mach number is large, ice shedding phenomenon will occur at the blade tip and the ice shape will change. In order to avoid the icing shedding problem, the condition of small Mach number is selected.

The grids are shown in Figure 2. The rotation period is divided into 120 time steps, and each time step requires 3000 iterations in flow field and droplet trajectory calculation, respectively. At least three complete cycles are used to establish the periodicity of flow field. The computing platform is Inter(R) Core(TM) i7-4790K CPU processor with a clock speed of 4.00 GHz and the memory is 16 GB. Each step of flow field and droplet trajectory calculation takes nearly 5 h. The ice accretion is calculated on the rotor surface, so it takes only a few minutes which can be ignored.

Droplet collection efficiency analysis

Figure 8 shows the droplet collection efficiency of rotor at different azimuth angles. Figure 9 shows the distribution of droplet collection efficiency on the $r / R = 0.9$ section at different azimuth angles. X indicates the chordwise distance of the airfoil, and the positive value represents the upper surface, the negative value represents the lower surface, and the 0 value is the leading edge. Figure 10 shows the variation trend of maximum droplet collection efficiency on the $r / R = 0.9$ section. The local velocity of rotor blade is influenced by the rotation speed and the advance ratio in forward flight, and the local velocity will vary at different azimuth angles. The change of local velocity will affect the magnitude and distribution of droplet collection efficiency. Therefore, the variations of local velocity have a great influence on the droplet collection efficiency. In the same way, the pitch angle also affects the results. It can be seen that the droplet collection efficiency is mainly distributed near the leading edge of the rotor blade, and the value and range of droplet collection efficiency vary with the azimuth angle.

Figure 8.

Distribution of droplet collection efficiency at different azimuth angles.

Figure 9.

Curve of droplet collection efficiency at different azimuth angles.

Figure 10.

Curve of maximum droplet collection efficiency.

The distribution of droplet collection efficiency changes with pitch angle. When azimuth angle is 90°, the pitch angle is the smallest. The droplet impingement zone is closer to the upper surface and the range of droplet collection efficiency is minimum. When azimuth angle is 270°, the pitch angle is the largest. The droplets mainly impact on the lower surface, and the range of droplet collection efficiency is the largest.

The maximum value of droplet collection efficiency is about 1, and the minimum is nearly 0.7. This is because the local velocity of blade leading edge varies with the azimuth angle in forward flight. When the azimuth angle is 90°, the local velocity of blade is the largest, and the droplet impact velocity also increases, so the droplet collection efficiency is the largest according to the calculation formula. When the azimuth angle is 270°, the leading edge has the smallest local velocity, so the droplet collection efficiency is the smallest. It can be seen that the variations of local velocity at different azimuth angles affect droplet collection efficiency. When the blade rotates in the windward half circle ( $ψ = 0 \circ - 180 \circ$ ), the droplet collection efficiency is larger, and when the blade is in the tailwind half circle ( $ψ = 180 \circ - 360 \circ$ ), the droplet collection efficiency is smaller. From the maximum droplet collection efficiency curve in Figure 10, it can be concluded that the droplet collection efficiency increases first and then decreases.

Icing characteristics analysis

Figure 11 shows the ice growth rate of rotor at different azimuth angles. Figure 12 shows the distribution of ice growth rate on the $r / R = 0.9$ section at different azimuth angles and the average ice growth rate in the rotation period. Figure 13 shows the variation trend of maximum ice growth rate on the $r / R = 0.9$ section.

Figure 11.

Distribution of ice growth rate at different azimuth angles.

Figure 12.

Curve of ice growth rate at different azimuth angles.

Figure 13.

Curve of maximum ice growth rate.

It can be seen that with the change of azimuth angle, the ice growth rate and icing range of blade change greatly. When the azimuth angle is 90°, the ice growth rate is the largest, but the icing range is the smallest. When the azimuth angle is 270°, the ice growth rate is the smallest, while the icing range is the largest. The variation trend of the maximum icing rate curve is similar to droplet collection efficiency, but the details are different. Icing is mainly affected by droplets impacted on the surface, so the icing characteristics are similar to the droplet impingement characteristics. However, icing is not only affected by the impacted droplets, but also by air flow. At different azimuth angles, the air flow field changes greatly, so the variation trend is slightly different.

The average ice growth rate in the rotation period is different from that at each step. The icing range is affected by the icing characteristics at each azimuth angle, so the icing range is the same as the maximum icing range in the rotation period. The red line in Figure 13 represents the maximum value of average ice growth rate, and it can be seen that the value is slightly smaller than the values at the azimuth angles of 0° and 180°.

The calculation results show that the variations of local velocity have a great influence on rotor icing. The droplet collection efficiency and ice growth rate of the advancing blade are larger, while that of the retreating blade are smaller. Therefore, the unsteady motion of blade must be taken into account in rotor icing calculation.

Finally, the three-dimensional ice shape is shown in Figure 14 and the ice shapes of different blade sections in forward flight are given in Figure 15. The icing occurs mainly at the leading edge of blade. As the spanwise distance increases, the icing thickness becomes thicker and the icing range is larger.

Figure 14.

Three-dimensional ice shape of rotor blade.

Figure 15.

Ice shapes at different rotor sections.

Influence of time step on rotor icing

The influence of time step on the ice shape is studied. The $N$ means that a complete rotation period is divided into $N$ time steps for calculation. In Table 2, the icing characteristics under different time steps are calculated and the influence of time step on the ice growth rate is given. It can be seen that when the time step $N$ is greater than 60, the calculation results change very little and the error is less than 5%. In this case, it is considered that the suitable time step is chosen. When $N$ is very small, the icing calculation has less samples and the error of result increases. The time step will affect the accuracy of rotor icing calculation in forward flight. It is necessary to choose a suitable time step to ensure the correctness of the icing calculation and computing process does not cost much time.

Table 2.

The effect of different time steps on icing results.

X coordinate	N = 120	N = 60	N = 30	N = 15	N = 8
0.0543	2.51E–09	2.51E–09	2.52E–09	1.86E–09	0
0.04922	2.02E–08	2.06E–08	2.10E–08	2.07E–08	9.57E–09
0.0341	9.47E–07	9.51E–07	9.57E–07	9.63E–07	9.52E–07
0.02426	4.77E–06	4.77E–06	4.80E–06	4.85E–06	4.68E–06
0.00967	1.88E–05	1.88E–05	1.87E–05	1.89E–05	1.85E–05
−4.60E–04	8.03E–05	8.04E–05	8.03E–05	7.99E–05	8.07E–05

Conclusion

In this article, a numerical simulation method of three-dimensional rotor icing in forward flight is presented. The unsteady calculation method is adopted: the overlapping grids are used and the blade grid is generated by pseudo three-dimensional method. The flow field is calculated using N-S equations in the ground coordinate system and the droplet trajectory is calculated by Eulerian model. Three-dimensional surface ice accretion model is applied to simulate the rotor icing, and the influence of centrifugal force on the ice shape is considered in the model. The ice growth rate at each azimuth angle in a rotation period is obtained, and the icing shape can be calculated by the average ice growth rate. Through icing simulation results, the droplet impingement characteristics, and icing characteristics of the rotor in forward flight are obtained. Further conclusions are drawn as follows:

The three-dimensional rotor icing in forward flight can be simulated using the numerical simulation method proposed in this article. Icing occurs mainly at the leading edge of blade, and the ice shapes of the different sections are not the same.

The forward flight and pitch motion of blade have a great influence on the icing process. The variations of local velocity and pitch angle will lead to the change of the magnitude and distribution of droplet collection efficiency and ice growth rate obviously. The values of droplet collection efficiency and ice growth rate are larger in the windward half circle ( $ψ = 0 \circ - 180 \circ$ ), but smaller in the tailwind half circle ( $ψ = 180 \circ - 360 \circ$ ).

The average ice growth rate in forward flight is significantly different from that at different azimuth angles. The unsteady method is needed in the calculation process, and the time step will affect the calculation accuracy.

Footnotes

Handling Editor: Mustafa Canakci

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the National Basic Research Program of China (“973” Project) under grant no. 2015CB755800 and the National Natural Science Foundation of China under grant no. 11402114.

ORCID iD

Zhengzhi Wang

References

Korkan

Dadone

Shaw

. Performance degradation of helicopter rotor in forward flight due to ice. J Aircraft 1985; 22: 713–718.

Coffman

Jr.

Helicopter rotor icing protection methods. J Am Helicopter Soc 1985; 32: 34–39.

Gent

Cansdale

. The development of mathematical modelling techniques for helicopter rotor icing. In: Proceedings of the 23rd aerospace sciences meeting (report no. AIAA-85-0336), Reno, NV, 14–17 January 1985. Reston, VA: American Institute of Aeronautics and Astronautics (AIAA).

Wright

. Validation results for LEWICE 3.0. In: Proceedings of the 43rd AIAA aerospace sciences meeting and exhibit (report no. AIAA-2005-1243), Reno, NV, 10–13 January 2005. Reston, VA: American Institute of Aeronautics and Astronautics (AIAA).

Narducci

Kreeger

. Analysis of a hovering rotor in icing conditions. In: Proceedings of the 66th annual forum and technology display, Phoenix, AZ, 11–13 May 2010, Phoenix, AZ: American Helicopter Society.

Bain

Cajigas

Sankar

et al . Prediction of rotor blade ice shedding using empirical methods. In: Proceedings of the AIAA atmospheric and space environments conference 2010 (report no. AIAA-2010-7985), 2010. Reston, VA: American Institute of Aeronautics and Astronautics (AIAA).

Cao

Zhong

. Tandem helicopter trim and flight characteristics in the icing condition. J Aircraft 2010; 47: 1559–1569.

Chen

Zhao

. Numerical simulations for ice accretion on rotors using new three-dimensional icing model. J Aircraft 2017; 54: 1428–1442.

Zhao

Chen

. New 3-D ice accretion method of hovering rotor including effects of centrifugal force. Aerosp Sci Technol 2016; 48: 122–130.

10.

Duque

EPN

Ahmad

. Helicopter rotor blade computation in unsteady flows using moving overset grids. J Aircraft 1996; 33: 54–60.

11.

Strawn

Duque

EPN

Ahmad

. Rotorcraft aeroacoustics computations with overset-grid CFD methods. J Am Helicopter Soc 1998; 44: 132–140.

12.

Kang

Kwon

. Unstructured mesh Navier-Stokes calculations of the flow field of a helicopter rotor in hover. J Am Helicopter Soc 2002; 47: 90–99.

13.

Wang

Zhao

Liang

et al . Numerical analysis on noise of rotor with unconventional blade tips based on CFD/Kirchhoff method. Chinese J Aeronaut 2013; 26: 572–582.

14.

Chen

Mccroskey

Obayashi

. Numerical solutions of forward-flight rotor flow using an upwind method. J Am Helicopter Soc 2012; 28: 374–380.

15.

Yang

. Efficient numerical techniques for simulating a rotorcraft flow field with overlapping grids. AIAA J 2015; 53: 1372–1383.

16.

Jameson

Schmidt

Turkel

. Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time stepping schemes. In: Proceedings of the 14th fluid and plasma dynamics conference (report no. AIAA 81-1259), Palo Alto, CA, 23–25 June 1981. Reston, VA: American Institute of Aeronautics and Astronautics (AIAA).

17.

Arnone

. Viscous analysis of three-dimensional rotor flow using a multigrid method. J Turbomach 1993; 116: 435–445.

18.

Shi

Zhao

Fan

et al . A new single-blade based hybrid CFD method for hovering and forward-flight rotor computation. Chinese J Aeronaut 2011; 24: 127–135.

19.

Spalart

Allmaras

. A one-eq18 turbulence model for aerodynamic flows. In: Proceedings of the 30th aerospace sciences meeting and exhibit, Reno, NV, 15–18 January 1992. Reston, VA: American Institute of Aeronautics and Astronautics (AIAA).

20.

Peng

Zhao

Zhu

. CFD calculations on the unsteady aerodynamic characteristics of a tilt-rotor in a conversion mode. Chinese J Aeronaut 2015; 28: 1593–1605.

21.

Bourgault

Boutanios

Habashi

et al . Three-dimensional Eulerian approach to droplet impingement simulation using FENSAP-ICE, part 1: model, algorithm, and validation. J Aircraft 2000; 37: 95–103.

22.

Honsek

Habashi

Aubé

. Eulerian modeling of in-flight icing due to supercooled large droplets. J Aircraft 2008; 45: 1290–1296.

23.

et al . Investigation of water impingement on a multi-element high-lift airfoil by Lagrangian and Eulerian approach. Propuls Power Res 2015; 4: 161–168.

24.

Thomas

Cassoni

Macarthur

. Aircraft anti-icing and de-icing techniques and modeling. J Aircraft 1996; 33: 841–854.

25.

Dong

Zheng

Zhu

et al . Calculation and analysis of runback water flow on anti-icing airfoil surface. J Aircraft 2016; 53: 1597–1605.

26.

Wang

Zhu

. Study of the effect of centrifugal force on rotor blade icing process. Int J Aerospace Eng 2017; 2017: 1–9.

27.

Caradonna

Laub

Tung

. An experimental investigation of the parallel blade-vortex interaction. Technical Manual NASA TM-86005, 1984. Moffett Field, CA: NASA Ames Research Center.

28.

Wirogo

Srirambhatla

. An Eulerian method to calculate the collection efficiency on two and three dimensional bodies. In: Proceedings of the 41st aerospace sciences meeting and exhibit (report no. AIAA-2003-1073), Reno, NV, 6–9 January 2003. Reston, VA: American Institute of Aeronautics and Astronautics (AIAA).

29.

Brouwers

Palacios

Smith

et al . The experimental investigation of a rotor hover icing model with shedding. In: Proceedings of the 66th annual forum proceedings—AHS international, Phoenix, AZ, 11–13 May 2010, pp.1157–1173. US: American Helicopter Society.

Numerical simulation for three-dimensional rotor icing in forward flight

Abstract

Keywords

Introduction

Numerical simulation method

Grid generation method

Air flow field calculation

Droplet trajectory calculation

Ice accretion calculation

Ice shape calculation

Method verification

Flow field verification

Droplet trajectory verification

Rotor icing verification

Analysis of rotor icing in forward flight

Droplet collection efficiency analysis

Icing characteristics analysis

Influence of time step on rotor icing

Conclusion

Footnotes

Declaration of conflicting interests

Funding

ORCID iD

References