Abstract
In view of the rotor icing problems, the influence of centrifugal force on rotor blade icing is investigated. A numerical simulation method of three-dimensional rotor blade icing is presented. Body-fitted grids around the rotor blade are generated using overlapping grid technology and rotor flow field characteristics are obtained by solving N-S equations. According to Eulerian two-phase flow, the droplet trajectories are calculated and droplet impingement characteristics are obtained. The mass and energy conservation equations of ice accretion model are established and a new calculation method of runback water mass based on shear stress and centrifugal force is proposed to simulate water flow and ice shape. The calculation results are compared with available experimental results in order to verify the correctness of the numerical simulation method. The influence of centrifugal force on rotor icing is calculated. The results show that the flow direction and distribution of liquid water on rotor surfaces change under the action of centrifugal force, which lead to the increasing of icing at the stagnation point and the decreasing of icing on both frozen limitations.
1. Introduction
The helicopters are most likely to encounter icing problems due to the lower flight height and slower flight speed. The helicopter rotor blade icing is easy to cause a flight accident because rotor is both the lifting surface and the operating surface. Therefore, the rotor icing problem needs to be studied for all-weather flight helicopters [1].
The main difference between rotor icing and fixed wing icing is that linear velocity of rotor blade varies along the spanwise direction, resulting in different droplet impingement range and mass on the rotor blade. At the same time, the liquid water movement process is affected by not only aerodynamic force and gravity but also the centrifugal force, eventually leading to the variation of ice shape.
From the 1940s, NASA, Sikorsky Aircraft, and other research institutions have carried out some researches on helicopter rotor icing problems including natural icing flight experiments, artificial icing flight experiments, and icing wind tunnel experiments. A large amount of experimental data has been accumulated, and many kinds of helicopters like Black Hawk have passed the airworthiness certification [2, 3]. With the development of computer technology, researchers have explored the use of numerical simulation technology to solve helicopter rotor icing problem. At present, LEWICE [4] software developed by NASA is the most popular in the research field of helicopter rotor icing numerical simulation. Its calculation method is as follows: the two-dimensional flow field results are obtained from the three-dimensional rotor flow field and then ice accretion is simulated based on fixed wing icing calculation method; finally the differential geometric method is used to expand the two-dimensional ice shape into three-dimensional. However, this method does not realize the three-dimensional rotor icing numerical simulation [5]. The other scholars have also carried out the related researches on the icing simulation of rotor blade airfoil, which is similar to the NASA research method, but its calculation accuracy still needs to be improved [6, 7].
In this paper, the effect of centrifugal force on helicopter rotor blades icing process is studied by numerical simulation method. A numerical simulation method is developed to simulate the blades’ ice accretion in rotating flow field. The Eulerian two-phase flow is used to calculate the droplet impingement characteristics. A method of water flow distribution based on shear force and centrifugal force is proposed to simulate three-dimensional blade ice shape. The calculation results and available experimental results are compared to verify the algorithms. Finally, the surface water film flow under centrifugal force is studied, and the effect of centrifugal force on the blade icing is analyzed.
2. Numerical Simulation Method
The icing numerical simulation includes the following modules: air flow field calculation, droplet trajectory calculation, and ice accretion calculation. Because of the small size of droplets, it is generally believed that the droplets’ movement is mainly affected by air flow, while the interference of droplets on the air flow field is neglected.
2.1. Air Flow Field Calculation
At present, overlapping grids [8, 9] are mainly used to calculate rotor flow field, which include the body grid and the background grid. Among them, three-dimensional grid with regular geometric shape is generally adopted in the background grid, so it is not difficult to generate the grid. The body grid is three-dimensional grid around the rotor blade. In this paper, the pseudo three-dimensional method is used to generate the body grid, which can reduce the computing time and resources and grid quality can meet the computational requirements. Pseudo three-dimensional grid generation is mainly divided into the following steps: first, a large number of cross sections are divided along the spanwise direction according to rotor blade shape; then two-dimensional grid generation method is used to generate grid of each section; the two-dimensional grid is connected to the three-dimensional grid based on the geometric topological relation; and finally grid of blade root and blade tips is generated. The overlapping grids used in this paper are shown in Figure 1. The rotor flow field in hover is axial symmetry, so the periodic boundary condition is adopted in the background grid. In order to accurately simulate the boundary layer, the first near-wall grid size of body grid is c and , where c is chord length.
In hover condition, the rotor blade is assumed to rotate around the shaft at a certain angular velocity and the blade flapping and lagging is not considered. The air velocity in the inertial coordinate system is the absolute velocity and in rotating coordinate system is the relative velocity . The air density , temperature , and air pressure are scalar, so the numerical value does not change in different coordinates. The Reynolds-averaged Navier–Stokes equations expressed by the absolute velocity in the rotating coordinate system can be presented as [10]where is position vector, is second-order unit tensor, is stress tensor, and is total energy per unit volume. The following relationship exists between the relative velocity, the absolute speed, and the rotational speed:It can be seen that when the rotation speed , the equations are equivalent to the N-S equations in inertial coordinate system.
In the calculation, the N-S equation in inertia coordinate system is solved in the background grid, and the N-S equation in rotating coordinate system is solved in the body grid. In order to make the data transmission between the two grids more convenient, the absolute velocity is used as the calculation parameter in the body and the background grid.
In the solver, the cell-centered finite volume scheme is used for spatial discretization. A second-order central difference scheme is employed in convective fluxes computation and the artificial dissipative terms developed by Jameson et al. [11] are added to prevent the appearance of oscillations. A five-stage Runge-Kutta scheme is used to solve the temporal discretization. In this work, the local time stepping and implicit residual smoothing are also applied to improve the calculation efficiency [12, 13]. One-equation Spalart–Allmaras (S–A) [14] turbulence model is adopted to simulate the viscous effect. S–A turbulence model is widely used in the rotor CFD calculation and the calculation accuracy can meet the requirement [8, 15].
2.2. Droplet Trajectory Calculation
For rotor icing, the Eulerian method is used to calculate droplet trajectory. The supercooled droplets are considered as a continuous phase in the Eulerian method and droplet volume fraction is used to represent the amount of liquid water within a control volume. The quality and velocity distribution of droplets in the calculation domain is obtained by solving the droplets governing equations. In the derivation of droplet trajectory equations, the following assumptions are made: the droplets are spherical with no deformation or breaking; the temperature of supercooled droplet is the same as the air temperature; therefore the heat and mass transfer between droplets and air are ignored and the physical properties of droplet are constant; there is no collision, coalescence, or splashing of the droplets; and air drag, gravity, and buoyancy are considered in the droplet movement [17].
Based on the air flow field results, the continuity and momentum equations for the droplets are listed as [18–20] where is droplet volume fraction, is water density, is droplet velocity, is air velocity, is acceleration of gravity, and is the inertial factor: is the droplet diameter; the aerodynamic coefficient; and is the drag function; in this paper it is represented as [17] is the drag coefficient; Re is the relative Reynolds number, which can be expressed asThe droplet collection coefficient is one of the most important parameters in the icing research, which is the ratio of the droplet quality collected and the maximum droplet quality. It can be defined aswhere is the far field droplet volume fraction, is the local droplet volume fraction, is the far field droplet velocity, is the local droplet velocity, and is the surface normal vector.
It can be seen that the continuity and momentum equations for the droplets are similar to the continuity and momentum equations for air, so the overlapping grids can also be adopted in droplet trajectory calculation and droplets governing equations are solved in the same manner as N-S equations. In this paper, the finite volume method is used to discretize the droplets governing equations, and the acceleration techniques such as local time stepping and implicit residual smoothing method are employed to speed up convergence.
2.3. Ice Accretion Calculation
The effect of centrifugal force is needed to be considered especially in the rotor ice accretion calculation. The following assumptions are made: the effect of centrifugal force on the liquid water movement is researched, but the effect on frozen ice is neglected. So ice shedding from the blade tip due to centrifugal force is not considered in this paper.
Based on the results of droplet collection coefficient, thermodynamic and fluid dynamic principles can be applied to calculate icing processes. In the ice accretion calculation, the droplets impinge on the rotor blade surface and gradually freeze. As shown in Figure 2, the blade surface is divided into a number of control volumes; the four adjacent volumes of each control volume are represented by W, S, E, and N, respectively. The mass and energy conservation equations of control volume are listed as [21, 22]In the mass conservation equation, is the mass flux of droplet collected:where is droplet collection coefficient, is the far field droplet velocity, and LWC is liquid water content.
is the mass flux of water flow into control volume:where is the mass flux of water from adjacent volume W flow into control volume, and so forth. The flow direction is consistent with Figure 2.
is the mass flux of water evaporation:where is convective heat transfer coefficient, is air density, is the water vapor gas constant, is air specific heat capacity at constant pressure, is the relative humidity, is surface temperature, is ambient temperature, and and is the saturated water vapor pressure of surface and far field, respectively.
is mass flux of water frozen; is the mass flux of water outflow from control volume. The amount of water outflow from control volume is equal to the amount of water flow into the next unit, and the calculation method is introduced in Section 2.4.
In the energy conservation equation, is the heat flux of the droplet collected including the internal energy and the kinetic energy; is the internal energy of water flow into control volume; is the internal energy of water outflow from control volume; is heat flux conducted by skin; because anti-icing is not considered in icing process, the surface is adiabatic, ; is convective heat transfer of water and air; and is the energy of water evaporation. The corresponding heat flux can be expressed as follows:where is isobaric specific heat capacity of ice, is isobaric specific heat capacity of water, is fusion latent heat of ice, and is latent heat of water evaporation. is phase equilibrium temperature of ice, is surface temperature, is static temperature, is surface recovery temperature, and is velocity of droplet impacted on surface.
is freezing coefficient which is developed by Messinger [23]. It is defined as the ratio of the mass flux of water frozen and the total mass flux of water in control volume:On the basis of the results of air flow field and droplet trajectories, the ice accretion equations can be solved. It is necessary to start the solution at the stagnation point where it is assumed that . Then all control volumes are solved according to the direction of water flow and can be obtained. The ice shape can be simulated by adding the ice thickness along the wall-normal direction, and the ice thickness is expressed aswhere is the icing time and is the area of control volume.
2.4. The Effect of Centrifugal Force
The movement of droplets impacted on the surface is influenced by the external force. In nonrotating state, the liquid water on surface is mainly affected by the shear force of air flow. Under the influence of shear force, the movement direction of liquid water is the same as the direction of shear force, that is, water flow from wing leading edge to trailing edge.
In the rotating state, the liquid water on surface is subjected to the combined action of shear force and centrifugal force. The magnitude and direction of shear force are the same as those in nonrotating state. The magnitude of centrifugal force is related to the mass of liquid water, the position of control volume, and the rotational speed, and the direction is always from blade root to tip, so the liquid water will flow along the spanwise direction under centrifugal force. Therefore, compared with nonrotating state, the magnitude and direction of the force are varied in the rotating state.
The shear force of air flow can be expressed asThe centrifugal force of each control volume can be calculated by the following formula:where is the shear stress for air flow, is the control volume area, is the local water film thickness, is the angular velocity, and is the spanwise length.
The resultant force in the rotating state can be expressed asFlow velocity can be expressed by the following formula:where is the control volume area, is the local water film thickness, and is the dynamic viscosity of water.
Compared to nonrotating state, the magnitude and direction of the resultant force change in the rotating state, leading to a variation of flow velocity and mass flux of liquid water outflow from control volume. According to the relevant research, the liquid water flow and mass distribution are mainly influenced by the shear stress and the centrifugal force [24]:The distribution principle of runback water in control volume is as follows:The runback water flow is decomposed into two directions of south-north direction (SN) and west-east (WE) direction. Take north-south direction as an example: when , there is no water flow; when , the liquid water flows from the control volume to the adjacent unit N; when , the liquid water flows from the control volume to the adjacent unit S.
and are the components of the resultant force in the WE and SN directions, respectively:Therefore, a distribution principle of runback water is proposed as follows:
3. Method Validation
In order to verify the correctness of the numerical simulation method used in this paper, rotor blade icing under three conditions is simulated and the calculation results are compared with the experimental results and the ARISP icing model in [16]. Blade airfoil is NACA0015, blade chord length is 0.15 m, and blade radius is 1.24 m. The calculation conditions are shown in Table 1.
Comparisons of ice shape under three conditions are shown in Figures 3–5. The icing temperature ranges from −5°C to −15°C, including three typical icing types, namely, the glaze ice, the rime ice, and the mixed ice. It can be seen that calculation results are in good agreement with the experimental results in three conditions. The predicted icing thickness is approximate to the experimental results, and the upper and lower surfaces frozen limitations are simulated better. And the calculation results are more close to the experimental results than the ARISP icing model. Therefore, under the three conditions of the glaze ice, the rime ice, and the mixed ice, the rotor blade icing can be simulated successfully using numerical simulation method in this paper, and the algorithm has better accuracy.
4. Results and Analysis
In order to further investigate the influence of centrifugal force on ice accretion during rotor icing process, numerical simulation method is used to calculate the liquid water distribution and ice shapes of rotor blades in the nonrotating and rotating state. The same flow field and droplet trajectory results are used in calculation, and the main difference between two states is that the effect of centrifugal force on ice accretion model is taken into account in the rotating state.
As analyzed in Section 2.4, the effect of centrifugal force on the rotor icing is mainly reflected in its effect on the liquid water flow. When the icing temperature is low, the droplet is frozen immediately after impacted on the surface, so the effect of centrifugal force is not obvious. When icing temperature is high, only part of water is frozen rapidly; another part will flow under the centrifugal force and aerodynamic force and phase change into ice. At this time, the influence of centrifugal force on ice shape is great. Therefore, condition1 is selected to investigate the influence of centrifugal force on the icing shape in this paper.
The calculation results of the blade ice shape under two conditions are shown in Figure 6. It can be seen that the ice shape calculated under the centrifugal force is more close to the experimental results. The centrifugal force in icing model has a great influence on the final ice shape: compared to the nonrotating state, the icing on leading edge increases and on both frozen limitations reduces. The ice thickness at the stagnation point increases from about 7 mm to 10 mm, which is 30% more than experimental results. This is because a part of liquid water changes its flow direction under the effect of centrifugal force. There is more liquid water flow along the spanwise direction, resulting in the water being frozen near stagnation point instead of the upper and lower frozen limitations.
In order to further analyze the reasons for the change of ice shape, the movement and distribution of liquid water on the surface are analyzed. Firstly, the coordinate system is established in Figure 7, in which the -axis is spanwise direction and the -axis is chordwise direction. Since the calculation airfoil NACA0015 is symmetrical, the control volumes on upper surface of airfoil are selected for analysis. Three positions are chosen along the spanwise direction, which are located in , and four positions are selected along the chordwise direction, which are located in m, 0.001 m, 0.002 m, and 0.003 m. There are 12 units in total and each unit area is 10−6 m2. The distribution and number of units are shown in Figure 7. Through the analysis of the force and flow of liquid water in 12 units, the influence of centrifugal force on the icing is obtained.
The liquid water movement is mainly influenced by the aerodynamic shear force and the centrifugal force, and the force of each unit is shown in Table 2.
In the nonrotating state, the surface liquid water movement is mainly affected by the aerodynamic shear force. It can be seen from Table 2 that the magnitude of aerodynamic shear force increases with chord length. Because the air flow through the leading edge is accelerated and the velocity increases gradually, the aerodynamic shear force increases and the direction is almost the same as the surface tangential direction. The aerodynamic shear force increases with spanwise position, which is due to the increase of span length and linear velocity.
The magnitude of centrifugal force is related to the mass of liquid water, the position of control volume, and the rotational speed. As we can see that, in the same spanwise position, the centrifugal force decreases with the increase of chord length, which indicates that the thickness of the liquid water is reduced, and the direction of centrifugal force is always from blade root to tip.
In the nonrotating state, the effect of centrifugal force makes the magnitude and direction of resultant force change. Aerodynamic shear force is relatively slight near the stagnation point ( m) and centrifugal force is larger due to the more liquid water mass. Therefore, the centrifugal force plays a dominant role in the resultant force and the direction of liquid water movement varies from chordwise to spanwise. However, aerodynamic shear force is relatively large near the impact limitation ( m) and the centrifugal force is less due to the decrease of surface liquid water mass, so the centrifugal force has little effect on resultant force and the direction of liquid water movement deflects slightly.
Liquid water distribution on blade surface is shown in Figure 8 and relevant data are given in Table 3. It can be seen that, in the nonrotating state, the thickness of liquid water increases first and then decreases with chordwise position. Although the water mass collected is larger in the leading edge, the surface liquid water is moving along the chordwise direction under the aerodynamic shear force, so the liquid water in the leading edge is relatively thin. Liquid water is gradually frozen in the moving process, and the typical glaze ice is formed.
(a) Nonrotating state
(b) Rotating state
The centrifugal force is considered in the rotating state. The flow direction of liquid water near leading edge changed from chordwise direction to spanwise, resulting in the liquid water thickness increase near the stagnation point and decrease near the frozen limitations. There is a direct relationship between the ice thickness and the liquid water distribution. Therefore, under the effect of centrifugal force, ice thickness increases near leading edge and decreases at the upper and lower limitations.
Above all, the effect of centrifugal force on icing is mainly reflected in the following: centrifugal force makes the magnitude of resultant force of liquid water increase and the direction change; liquid water is trended to flow along spanwise direction, which lead to the thickness of liquid water and ice accretion increase near the leading edge and decrease at the upper and lower frozen limitations.
5. Conclusions
In this paper, numerical simulation method is used to study the influence of centrifugal force on rotor blade icing. A three-dimensional numerical simulation method for rotor icing is established and a runback distribution method based on the shear force and centrifugal force is proposed. The three-dimensional ice shapes of rotor blades are simulated in different conditions and the calculation results are compared with the experimental results in different conditions include glaze ice, rime ice, and mixed ice. The influence of centrifugal force on icing process is analyzed qualitatively. The following conclusions are obtained:(1)Three-dimensional rotor blade icing can be simulated using the numerical simulation method presented in this paper. The calculated results are in good agreement with the experimental results, which proves that the method is effective.(2)The effect of centrifugal force on the liquid water is obtained through numerical calculation method. It can be seen that due to centrifugal force, the magnitude of resultant force is larger and direction is changed, resulting in liquid water distribution change.(3)Centrifugal force plays an important role in the rotor icing process. Centrifugal force makes the ice shape of rotor blade change. Ice thickness increases near the stagnation point and decreases at the upper and lower frozen limitations.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Basic Research Program of China (“973” Project) (Grant no. 2015CB755800), the National Natural Science Foundation of China (Grant no. 11402114), and Jiangsu Province Scientific Research Innovation Project (Grant no. KYLX-0314).