Abstract
The main objective of this study is to estimate the dynamic loads acting over a glaze-iced airfoil. This work studies the performance of unsteady Reynolds-averaged Navier-Stokes (URANS) simulations in predicting the oscillations over an iced airfoil. The structure and size of time-averaged vortices are compared to measurements. Furthermore, the accuracy of a two-equation eddy viscosity turbulence model, the shear stress transport (SST) model, is investigated in the case of the dynamic load analysis over a glaze-iced airfoil. The computational fluid dynamic analysis was conducted to investigate the effect of critical ice accretions on a 0.610 m chord NACA 0011 airfoil. Leading edge glaze ice accretion was simulated with flat plates (spoiler-ice) extending along the span of the blade. Aerodynamic performance coefficients and pressure profiles were calculated and validated for the Reynolds number of 1.83 × 106. Furthermore, turbulent separation bubbles were studied. The numerical results confirm both time-dependent phenomena observed in previous similar measurements: (1) low-frequency mode, with a Strouhal number –0.02, and (2) higher frequency mode with a Strouhal number –0.69. The higher frequency motion has the same characteristics as the shedding mode and the lower frequency motion has the flapping mode characteristics.
1. Introduction
Ice accretion on aerodynamic surfaces leads to significant deterioration of the blade performance and operation [1]. Heavy icing disrupts the continual power generation from the wind turbines and the disruption may be prolonged in severe cold conditions. The database published by “Statistics Sweden” showed 161,523 hours of total downtime for the period 1998–2003 where 7% of the downtime was related to icing condition and resulted in over 5% production loss in the country [2].
In 1958, the effects of ice formation on different airfoil sections were studied by NACA [3]. They measured lift, drag, and pitching moment coefficients of an NACA 65A004 airfoil section. The icing on the airfoil caused a rapid increase of the drag coefficient, a drop of the lift coefficient, and change in the pitching moment coefficient. In 2012, Villalpando et al. [4] conducted numerical simulations over a two-dimensional ice-accreted NACA 63-415 airfoil at various angles of attack. They validated the load coefficients with experimental data at one angle of attack and extracted the modified pressure distribution due to the ice accumulation. The ice accumulated over the airfoil modified the pressure distribution affecting considerably the aerodynamic performances. Further, in 2014 [5] experimental and numerical investigations demonstrated that with increasing angle of attack the degradation of the instantaneous lift coefficient increases with a linear process.
Ice accretion is a phenomenon where super-cooled water droplets impinge and accrete on a body. Icing occurs on the leading edge of wind turbine blades [6]. The ice likely to form on wind turbine blades is of two kinds [7]. The first is called “rime ice,” which is generated at very low temperature, below –10°C. In rime ice conditions, the droplets in the air freeze instantly at the impingement point. The second is called “glaze ice,” which is generated in the temperature range −10 to 0°C. On glaze ice conditions, droplets gradually freeze while moving along the body, so-called runback [8].
Glaze ice accretion is often characterized by the presence of large protuberances, commonly known as glaze horns (Figure 1), which can cause flow separation downstream of the horns [9]. On iced airfoils, the boundary layer separates near the top of the horn, due to the pressure gradient produced by the large discontinuity in the surface geometry. This is a reason why CFD prediction is more challenging in the case of airfoil/blade with glaze ice than with rime ice [10].
In 1986 [11] experiments were conducted to study the aerodynamic characteristics of the NACA 0012 airfoil under glaze ice accretion. Lift and drag penalties due to the ice shape were found and the surface pressure showed the presence of large separation bubbles. A large flow separation region was observed and correlated to the pressure measurements. The iced cases were also analyzed numerically at angles of attack below stall [12]. In 2018, the ice accretion transient phenomenon and its effect on turbine performance was studied by coupling 2D steady state CFD, with a blade momentum and an ice accretion code [13].
As the specification of the horn is determinative in the subsequent vortices structure, a parametrization of these horns is desired. Although the exact size and shape of the ice formations are found to be complex functions of both the operating and icing conditions, a horn shape can be characterized by its height, the angle it makes with respect to the chord line, and its location [14]. In the mentioned studies ([11, 12]), the simulated ice shapes were constructed to approximately duplicate an actual measured ice accretion. The ice was accreted in an Icing Research Tunnel on a NACA 0012 airfoil. In an experimental study [15], glaze ice shapes were simulated by means of spoilers attached on the upper surface of an airfoil near the leading edge as a single-horn glaze ice. The spoilers modeled in this work were sized to simulate 22.5-minute glaze ice accretions for the NACA 0011 airfoil. Leading edge glaze ice accretions characteristically may consist of an upper and lower surface horns, called double-horn ice. In 2000, experimental data were provided for the NACA 0011 airfoil with simulated glaze ice shapes on both upper and lower surfaces [9]. The ice accretions can also be simulated with wooden forward-facing quarter-round shapes. Lee and Bragg [16, 17] have tested some geometries consisting of backward facing quarter-round, half-round, and forward-facing ramp. In [18], the simulated glaze horn-type ice accretions were determined from averaging geometry data from a set of actual ice accretions collected in a test at an Icing Research Tunnel. Single-horn simulations were used for that research by means of a 3 by 3 matrix of ice shape size and radius which was designed to parametrically vary these parameters. In [19] a single-horn glaze ice was placed on the upper surface of the airfoil which was a combination of a semicircle and a rectangle. In the previous study [20], steady state simulations were performed to investigate the mean flow characteristics. The aerodynamic performance coefficients and pressure profile were determined with CFD and compared with the available measurements in [9]. The steady state simulations were not found to be reliable to calculate the loads in the case of iced airfoil. Time-dependent simulations were conducted to determine the vortices structure. It was concluded that the thickness of the leading edge glaze ice should be taken as an effective parameter of an ice shape. Although it does not change the formed vortices on the airfoil downstream of the spoiler tip, it can affect the acting forces on the spoilers. As a matter of fact, the main part of the total lift force on the airfoil was found to be on the spoilers where the ice shape may significantly affect the lift and drag.
The leading edge ice accretion not only may be detrimental to aerodynamic performance but also can be a concern in terms of large unsteady loads associated with the flow separation [21]. Generally, the steady state effects of the separation bubbles on the airfoil performance are characterized by large increases in drag, reductions in lift, and changes in airfoil pitching moment characteristics. The separation occurs at the tip of the horn ice shape (Figure 1) and immediately a shear layer forms, which separates the recirculation region from the freestream flow. The shear layer begins to roll as it moves downstream, and vortices within the shear layer merge, forming larger vortical structures.
Understanding the behavior of the separation bubble is critical to understand the effect of ice accretion on airfoil aerodynamics. Laminar separation bubbles on airfoil have been widely studied [22–26]. The bubble forms when a laminar boundary layer encounters an adverse pressure gradient of sufficient strength to cause separation. The separated flow may be divided into two main regions: The free shear layer and the recirculation bubble. These two regions may then be further subdivided into parts, upstream and downstream of the transition point. After transition, the magnitude of the reverse flow increases and a vortex type flow is seen in the bubble. Before the transition, the reverse flow is very slow and this area is sometimes referred to as a dead-air region.
Flows at high Reynolds number with separation and reattachment have long been a subject of many studies [27]. In 1983, the structure of a turbulent separation bubble was studied by Kiya and Sasaki [28]. In their measurements, they observed a large-scale unsteadiness accompanied by an enlargement and shrinkage of the bubble and a flapping motion of the shear layer near the separation line. At a turbulent separation bubble, the flow is characterized by two separate time-dependent phenomena: flapping and shedding [29]. One mode is associated with a global breathing motion of the separation bubble, described as “flapping motion” in the literature. The other mode is associated with the roll-up of spanwise vortices in the shear layer above the recirculating region and their shedding downstream of the separated zone. The regular mode includes the presence of vortical motion in the separated shear layer and vortex shedding from the separation bubble. The source of the regular mode has been attributed to the Kelvin-Helmholtz instability, where the difference between the velocity within the recirculation region of the separation bubble and the external flow causes a roll-up and shedding of vortices within the shear layer [30]. The shear layer flapping is characterized by low frequencies at certain locations in the separation bubble downstream of the separation onset. The oscillations tend to operate across different frequency scales compared with those of the regular mode [30].
As the shear layer vortices are formed and shed, the height and length of the separation bubble tend to change as a function of time [11]. The boundary layer events produce variations in lift, drag, and moment coefficients [31]. It has been found that the frequency of the oscillating flow can be nondimensioned with the Strouhal number St based on the momentum thickness [31], airfoil projected height [21, 26, 30], or separation bubble length [30], where is the flow oscillation frequency, is the separation bubble length, is the airfoil chord, θ is the boundary layer momentum thickness, α is the angle of incidence, and is the freestream mean velocity. Most results from the literature identify the effects of shear layer flapping as occurring at a Strouhal number on the order of or . The regular mode reported in the literature consistently corresponds to a Strouhal number in the range of 0.5 to 0.8 [30].
In 1987, Rumsey [32] used a numerical method to predict unsteady flow over different airfoil geometries at high angles of attack. Using a compressible, two-dimensional, Navier-Stokes code, Rumsey computed the flow over a NACA 0012 airfoil, without any imposed perturbation, at a Reynolds number and a Mach number . A low-frequency oscillation in the flow was encountered, if a turbulent boundary layer near the leading edge was assumed. For laminar flow at , was found to be independent of the angle of attack over 20 degrees at a constant value of about 0,155.
In 1989, Zaman et al. [33] studied the low-frequency oscillations in the flow over NACA 0012 before stall. A turbulent boundary layer was resolved in a two-dimensional Navier-Stokes code. Details of the flow field and unsteady forces compared reasonably well with the experimental data. This study was explored experimentally as well as computationally for NACA 0012 airfoil with a “glaze ice accretion” at the leading edge [21]. With a Navier-Stokes computation, “limit-cycle” oscillations in the flow and in the aerodynamic forces were observed at low Strouhal number. They found that the occurrence of the oscillation depended on the turbulence model. With respect to the computations, questions remain in the application of turbulence models to separated flows. Nevertheless, they concluded that the essence of the phenomenon can be captured computationally with certain combinations of the turbulence model, Reynolds number, and airfoil shape [21].
The separation bubble on an airfoil at low Reynolds number behind a simulated leading edge glaze ice accretion was studied experimentally in 1992 [26]. Time-dependent measurements of the flow field were performed for a laminar separation bubble. In 2002 [34], unsteady pressure measurements were performed on NACA 0012 airfoil with two- and three-dimensional leading edge glaze ice accretions. The mean and fluctuating lift coefficients at different angles of attack were presented. Gurbacki noticed the formation of additional vortices in the separated shear layer due to the jaggedness of the three-dimensional ice shapes. The unsteady content of the iced airfoil flow field was further analyzed [35]. The iced airfoil performance and distributed surface pressure were found to be like the unsteady 2D separation bubble in simple geometries. The NACA 0012 airfoil was tested at a high Reynolds number in 2013 [36]. In addition to the clean configuration, the airfoil model was also tested with a set of boundary layer trips, a two-dimensional extrusion of a horn ice shape casting, and an array of simulated icing configurations created using simple geometries. The resulting values of Strouhal number exhibited a dependence on the airfoil angle of attack and corresponded to a range that was consistent with the Strouhal number values reported in prior studies of the low-frequency mode in the literature.
With the above background, the various aspects of the separated flow over a glaze ice have been explored experimentally from 1986 until recently, while the capacities of the numerical computational methods in this area need more clarification. This work studies the performance of URANS simulations in predicting the oscillations over an iced airfoil. It is of interest to see how accurate URANS models estimate the structure and the size of the formed time-averaged vortices in comparison to the measured ones. Furthermore, the accuracy of a two-equation eddy viscosity turbulence model, the SST model, coupled with a URANS model to determine the dynamic load over a glaze-iced airfoil needs to be investigated.
The work presented in this paper is a continuation of a previous work [20], which aimed to determine the aerodynamic performances of wind turbines in icing condition. The main objective of this study is to develop a numerical model to observe and quantify the effect of the unsteady flow over the modeled iced airfoil in the presence of glaze ice. It will help to understand the mechanism by which the dynamic loads are initiated and sustained and their magnitude. The numerical simulations capture the separated shear layer and the structure of large-scale vortices in a turbulent separation bubble.
To this purpose, time-dependent simulations of the spoiler-ice test are performed at different angle of attack [9]. The results are validated with experimental data. Then, the load fluctuations are analyzed at different angle of attack. The frequency of the load cycle is studied as well as the extreme values of the loads acting on the airfoil.
2. Materials and Methods
The geometry and boundary conditions are considered from the literature [9, 15]. The experiments were conducted in a single-return closed circuit wind tunnel for a chord Reynolds numbers of 1.83 × 106. The test case was a 0.610 m NACA 0011 airfoil with the spoilers on the upper and lower surfaces.
The experimental mean pressure distribution over the blade surfaces is available for different angles of attack and different spoiler installed angles as well as the lift coefficient, drag coefficient, and pitch moment coefficient.
The fluid problem is solved with the finite volume technique using the CFD code ANSYS CFX 15.0 solver in which the set of equations are the unsteady Navier-Stokes equations in their conservation form [37]. Transient simulations were conducted on the 2D geometry of the NACA 0011 previously described.
High-resolution advection scheme was selected for the spatial discretization and second-order backward Euler scheme was applied for the temporal discretization of the equations. The time step size is set to seconds based on a time step size analysis prescribed in [20].
SST model with automatic wall function was activated to model the turbulent flow. The - SST model appears to be an accurate turbulence model for boundary layer detachment prediction [31, 38]. Symmetric boundaries in a 2D model, as used for these computations, limit the application of more detailed model such as Detached Eddy Simulation model or Scale-Adaptive Simulation [39].
The computational domain has the same dimensions as the test section of the wind tunnel, 2.13 m high, 3.05 m wide, and 3.66 m long. The spoiler angle is set to −40° and 0° for the upper and lower one, respectively (Figure 2). To approximate the horn heights of 22.5 min glaze ice accretions, the spoilers are 3.81 cm. Because the blade geometry and the spoilers have a constant cross-test section, symmetry boundary conditions are used at both sides of the created 2D model. The spoiler walls boundary conditions are defined as no-slip walls as well as the top and bottom walls of the wind tunnel. Flow is assumed to have a uniform velocity at the inlet and is allowed to move backward at the outlet.
The convergence criterion was set to a root-mean-squared value of 10−6. The simulation continued until a periodic variation was achieved for the drag, lift, pitching moment, and surface pressures. Due to the lift force fluctuations with time, a large number of iterations for the convergence of the mean value are needed. Computations were performed in order to assess the ability of the present method to accurately predict steady and unsteady airfoil behavior both below and above maximum lift conditions.
The surface pressure is computed on the upper and lower surfaces of the airfoil, from the leading edge to the trailing edge. In one cycle of the load oscillations, the pressure distribution is recorded, and the averaged value is used as the mean pressure distribution. The dynamic pressure ) is used to nondimension the mean pressure and to define the pressure coefficient , where is the static pressure.
For each time step, lift and drag forces are calculated as the loads over the iced airfoil, including the spoilers. The lift and drag coefficients ( and ) are defined by dividing the respective force by the reference area and . The reference area is calculated from the airfoil chord length and the span thickness, which in a 2D case is the thickness of the layer cell that is modeled. The pitching moment is calculated with respect to the quarter chord location as reported in the experimental tests. The moment is normalized using , , and .
In the calculation of the Strouhal number, the spoiler-iced airfoil projected height at each angle of attack is considered as “” in . Because the separation bubbles cover almost the entire blade surface from the spoilers to the trailing edge, the chord length is also considered as a rough magnitude of the bubble length “” in .
2.1. Mesh Analysis
A multiblocking mesh consisting of hexahedral elements is generated. Grid topology of C-type was used to generate high quality structured hexahedral grid using the code ICEM CFD. The value of corresponding to the first grid point above the walls is set below one. The mesh is denser on the spoilers and the airfoil surface to resolve the boundary layer formed on these parts.
Before the simulations were performed at different angles of attack, mesh scaling tests and mesh performance tests were carried out. After the mesh scaling test, a total of 1.7 million hexahedral mesh elements were created. A preliminary mesh containing 0.41 million hexahedral elements was generated and refined around the airfoil and spoilers to maintain . Both sides of each spoiler need a different level of refinements as the dimensionless distance depends on the wall shear stress as well as the size of the first cell normal to the wall. The wall shear stress value is found different at each side of the spoilers as the flow is different on each side.
The mesh analysis is performed keeping the initial condition and the number of processors used. However, these parameters have shown some effect on the numerical solution. The mesh resolution is studied by resolving the mesh uniformly in all directions. Load coefficients have been considered as the key parameters. Four grids were used consisting of 0.41, 0.93, 1.7, and 3.9 million elements.
The simulation results are plotted in Figure 3. The accumulated value of the lift converges after a large number of iterations. The drag value was not found sensitive to the mesh density, not plotted here. The numerical uncertainty due to the mesh density is calculated based on three mesh densities (Table 1) as suggested in [40]. Three different grid densities (fine, ; medium, ; and coarse, ) were used for the scaling test.
Table 1 shows the computed parameters based on the procedure described in [40] to determine the GCI. 2% difference between the medium and fine grid results was observed. The converged solution obtained using the medium grid (1.7 M) was used in the subsequent simulations at different angle of attack.
3. Results and Discussion
3.1. Mean Value Study
Mean pressure distributions are shown in Figure 4 and compared to the experimental results of Papadakis et al. [9]. At the leading edge , the pressure coefficient is greater than one (1.06) for the simulations in all angles of attack. The 6% extra amount of energy comes from the content of static pressure at the inlet, in addition to the dynamic pressure. Assuming a zero-pressure flow at the exit of the experimental test section, a nonzero pressure (~150 pa) is necessary to overcome the viscous losses with a uniform velocity of 45 m/sec.
(a) °
(b) °
(c) °
In a clean airfoil at °, just after the stagnation point the flow speeds up on both upper and lower surfaces. The maximum pressure is reached at the leading edge as the stagnation point, and then the pressure switches to the suction passing through the nose of the airfoil. The term “suction” is used to indicate a pressure lower than the reference pressure. Regarding the airfoil camber, the maximum suction occurs around in which the flow has the highest velocity. Then, the velocity decreases, restoring the pressure at 85% of the chord length, and leaves the trailing edge with a positive pressure.
When it comes to the iced airfoil, the flow acceleration is different at the nose of the airfoil. Since the pressure distribution does not follow the pattern of the pressure distribution on the clean airfoil, the flow is not following the curvature of the airfoil; that is, it is separated. Comparing to the clean case, more suction emerges on both sides (upper and lower surfaces). It also confirms that the flow is following a path with a bigger curvature than the airfoil nose, which is a vortex perimeter. Based on the experiments, the flow speeding continues up to and 0.1 for upper and lower surfaces, respectively. The pressure distributions are not symmetrical anymore and the upper surface is under higher suction. Then the suction magnitude of the flow starts to remain constant. It shows that the flow has just passed the corner of the oval shape vortex and is going to traverse the above part of the vortex. As long as the pressure coefficient is constant, the thickness of the oval vortex is unchanged.
Although the flow reaches almost the same amount of suction as the experiments, the simulated flow separates immediately at , where the flow encounters the spoilers. It can be concluded that the simulation captures accurately the height of the vortex but the maximum thickness of the vortex is obtained closer to the leading edge. As the pressure remains constant, the vortex shape does not have any curvature and the flow is passing over a straight path. At around 35% of the chord length, the flow starts to decelerate on the airfoil upper surface. This shows that the flow is descending the right-side curvature of the oval vortex which is getting thinner. In the simulations, this phenomenon occurs at around 45% of the chord which means the width of the vortex is overpredicted. Experiment shows that the flow is reattached at around as it is following the same slope of the pressure distribution on the clean airfoil. In the simulated model, the flow does not reattach completely before getting to the trailing edge. Close to the trailing edge, the flow follows a similar path to the airfoil that shows the separated region has become narrow though not reattached yet.
On the lower surface, the flow reattaches at around , while in the simulation a long bubble extends to the airfoil trailing edge with a constant pressure coefficient. The simulated vortex is like a wide oval with a uniform height from leading edge up to 80% of the chord length. The uniform shape of the vortex can relate to the zero angle of the spoiler on the lower surface , while on the upper surface the vortex tends to follow the bend of the spoiler to the left and so the pressure coefficient does not remain constant. In the last 3-4% of the chord length , a vortex emerges just close to the trailing edge on both sides of the airfoil surface. That is the vortex shedding which transmits downstream and will be discussed later in this paper (see Figure 10).
At °, more suction magnitude is seen on the lower surface than the upper surface. It is in opposite to the flow at °. For the clean airfoil at a positive angle of attack, the lower surface is windward. So, the leeward flow on the upper surface has more suction. But in the case of iced airfoil, the separated flow reaches more suction on the lower surface at the separated region. Therefore, a negative lift is produced at 58% of the chord length based on the experiments (from to 0.65). The lift is then positive until the last 4% where the trailing edge vortex is formed. The simulated vortex on the lower surface starts further than in the experiments but has the same form when it continues. The vortex on the upper surface is reattached at while it extends to the trailing edge in the simulated flow.
At °, there is no suction on the lower surface of the clean airfoil, while a high suction is seen on the leeward upper surface. That is the source of a large positive lift at this angle of attack. When it comes to the iced airfoil, there is a negative lift from the leading edge up to and then the lift direction changes to positive, which is further discussed later.
At °, the flow on the lower surface does not reach the same suction as in the experiments. It means that a smaller vortex is modeled in the simulation. Although it is formed further down and its thickness is smaller than the experiments, the form and the curvature of the oval vortex are modeled similarly. Regarding both experiments and simulation results, the flow will not reattach on the lower surface and it extends to the trailing edge. On the upper surface, there is a wide uniform bubble that extends almost all over the airfoil length up to the trailing edge vortex. More suction is seen in the simulated vortex which means the vortex size is overestimated. Furthermore, it turns at and then the trailing edge vortex forms at which is also much bigger than the experimental results.
The load coefficients are shown in Figure 5 regarding the mean value of the loads. It will be discussed in the following section.
(a) Lift coefficient
(b) Drag coefficient
(c) Pitching moment coefficient
3.1.1. Lift Coefficient
Lift performance for the clean and iced airfoil is presented in Figure 5(a). Although in a clean airfoil, a higher angle of attack results in more lift; a different trend is seen in the case of glaze-iced airfoil.
For the clean airfoil, the angle of stall is 15° with corresponding maximum lift coefficients of 1.4 for Reynolds number of 1.83 million. The lift coefficient remains linear with the angle of attack up to approximately 10° [9]. The 3.81 cm spoiler-ice shapes at 2% chord on the upper and lower surfaces result in a considerable change in the lift coefficient. These changes include lift sign reversal at low angles of attack while large reductions in the lift and a reduction in the lift slope are observed at higher angles of attack.
The lift sign reversal can be explained by the pressure distribution presented in Figure 4. At °, regarding the experiments, in most of the chordwise parts there is more suction on the upper surface than the lower surface which results in a positive lift (Figure 4(a)). From to the trailing edge, although there is more suction on the lower surface, the difference between the lower and upper surfaces is very small. So, there would be a positive lift on the entire airfoil (Figure 5(a)). In the simulated flow, the bubble on the lower surface has an oval shape with almost uniform height since the suction magnitude is almost constant , while in the experiments the oval-shaped bubble seems to be thinner near the extremities. It leads to a higher negative lift on the later 35% of the chord length compared to the experiments. In summary, the lift on the airfoil is around zero based on the simulations, though still positive (Figure 5(a)).
At °, there is a negative lift in the first 57% of the chord length based on the experiments (Figure 4(b)). Then, the reattachment process on the lower surface decreases the suction magnitude and a positive lift is obtained, while the negative part is dominant and the net lift on the airfoil is negative (Figure 5(a)). In the simulated flow, the bubble on the upper surface is more flat compared to the experiments and it keeps the constant suction of . Thus, more positive lift is obtained over the later 40% of the airfoil length that compensates for the negative lift over the rest of the airfoil and the net lift is around zero. The curvature of the clean airfoil is designed to generate a very low suction on the lower windward side at ° (Figure 4(b)), which is increased due to the accreted ice. Also, the slope of suction decrease is different on the upper leeward side at the iced airfoil. It shows that the flow does not follow the curvature of the airfoil; that is, it is detached. It results in more lift loss comparing to zero angle of attack. At °, the high pressure on the windward side of the clean airfoil (Figure 4(c)) results in a big lift magnitude (Figure 5(a)).
At the iced airfoil, flow separation leads to lift sign reversal on the first 47% of the airfoil (Figure 4(c)). Then the flow reattachment processes on the lower side decrease the suction. Since the upper surface is still on a constant suction, a positive lift is obtained over this part of the airfoil which dominates the negative lift, and the net lift is positive. More lift in the simulated flow (Figure 5(a)) comes from the fact that the vortex on the upper surface is overestimated leading to a higher level of suction () compared to the experiments () (Figure 4(c)). It results in a larger difference between the upper and lower surface suctions, that is, more positive lift.
It seems that although the lift is increasing with increasing the angle of attack above 4° in both clean and iced airfoil, lift loss is too much due to the spoilers.
3.1.2. Drag Coefficient
For the clean 0.61 m NACA 0011, the lowest drag coefficient was in the range of 0.008–0.009 near °. Near stall, the drag coefficient reached a value of about 0.03 as shown in Figure 5(b). For the iced airfoil, the increase in drag was in the range of 1000% to 2000% with respect to the clean airfoil for angles of attack between 0 and 14 degrees [9]. In the simulated flow, the drag is even larger due to the overestimation of the bubble sizes as described.
In both the iced airfoil and the clean airfoil, the drag increases when the angle of attack becomes larger than 4°. However, the slope is higher in the situation of the iced airfoil. Therefore, it is not justified to have an iced airfoil operating at a higher angle of attack because the gain in the lift implies a rapidly increasing drag. In the case of the clean airfoil, a higher angle of attack (up to the stall) would be recommended since it would provide more lift without extra drag.
The drag coefficient is almost the same at ° and °, while it increases at °. Looking at Figure 4, it can be seen that in all of the three cases there is a bubble extending all through the airfoil surface (with a constant suction) and another bubble that becomes thinner in the reattachment process. Considering the curvatures of the pressure distributions, the first bubble is similar at ° and °, while it is less elongated and reattaches faster in the case of °.
3.1.3. Pitching Moment Coefficient
Clean airfoil pitching moment coefficients at 25% of the chord length location are presented in Figure 5(c). Positive pitching moment was observed for angles of attack in the range 0 to 15 degrees, indicating that lift force was acting ahead of the 25% chord point [9]. The variations are very small when the angle of attack is increasing. The effects of the spoiler-ices on the pitching moment characteristics of the airfoil are presented in Figure 5(c). A large deviation from the clean airfoil pitching moment is observed, including a sign reversal throughout the range of angle of attack.
The gap between the simulated flow and the experiments comes from the different positions of the vortices as previously described. Based on the pressure distributions, the different positions of the formation and reattachments of the vortices lead to a different force distribution with respect to the quarter chord point. This is the reason for the gap between the experiments and the simulations, while both the experiments and simulation results show that a higher angle of attack leads to a large pitching moment coefficient which is not desired.
From all the investigated cases, the simulation performed for the largest angle of attack, °, showed the least agreement to the experimental data. Figure 7 helps to understand the flow behavior at this angle of attack. The turbulence kinetic energy is shown for ° and ° as it is a representative of the turbulence intensity or turbulence level . Both plots are extracted at a time step in which a high positive lift occurs (6.7 for ° and 5.9 for °). The turbulence level looks higher at a higher angle of attack. It is seen that a large turbulent wake arises at the trailing edge which is originated from the upper surface of the blade at °. This wake is not damped at the modeled downstream domain.
The pressure distribution also confirms the presence of a large turbulent wake (Figure 4(c)). The trailing edge vortex shedding forms at which is farther from the trailing edge compared to the cases ( at ° and at °). The vortex induces a suction magnitude of that is twice the suction magnitude of the trailing edge flow at °, for instance. Therefore, the shedding vortex is wider and thicker for this case which leads to a higher turbulence intensity at the trailing edge and consequently downstream.
The presented computational versus experimental mean flow analysis enables validating the current simulation. The following results are from computations only and are assumed to reproduce flow physics correctly. The general characteristics are validated towards the features of similar iced airfoils reported in the literatures.
3.2. Time-Dependent Study
Usually, flow field studies of separation bubbles are focused on the time-averaged characteristics. However, the bubble flow fields are known to have strong unsteady characteristics that also play a role in the aerodynamic characteristics [14]. These unsteady features are now discussed.
The streamlines around the iced airfoil illustrate that the vortices form and move over the airfoil surface. Consequently, the pressure on each point of the airfoil surface changes with time as the vortices are formed and convicted. It causes the lift to oscillate from positive to negative repeatedly, leading to a dynamic loading of the airfoil.
The instantaneous pressure values are shown in Figure 7 for some random points on the lower and upper surfaces of the airfoil upstream and downstream of the spoilers. The plots represent the pressure on these points on the airfoil at a zero angle of attack during 0.07 s of the simulation time, consisting of around 2000 iterations (the simulation time is normalized by the time step value, i.e., 3.5 × 10−5). The frequency of the pressure is almost the same for all the points. This frequency is illustrative of the roll-up of spanwise vortices in the shear layer above the recirculating region and their shedding downstream of the separated zone. Closer to the trailing edge the amplitude of the pressure increases. This can be related to the vortices shedding downstream of the trailing edge.
Upstream of the spoilers (points 1, 2), there is a high positive pressure which is almost constant during the time. As it will be later shown in Figure 11, these points are located in the upstream quasi-stagnant vortices; thus the pressure is high and the time variations of the pressure are very small. Downstream of both spoilers, the pressure is negative as there is a suction underneath the formed vortices.
It is observed that for each pair of points with the same chordwise location, that is, 3-4 and 2–5, the pressure oscillations are out of phase. This indicates that the vortices are forming and leaving the airfoil surface in a periodic manner.
Because of the pressure fluctuations on the airfoil, the loads fluctuate. The integrated load time history is shown in Figure 8 for the simulated cases. The drag variation is smaller than the lift variation. It indicates that the vortices distribution is shifting more frequently between the upper and lower surfaces, affecting mostly the vertical loads. However, both drag and lift oscillate at a similar frequency.
At °, a bimodal frequency pattern is seen in the drag force, with the maximum drag varying cyclically between two values (at two directions), while there is a quasi-sinusoidal trend of the variations for the other cases; that is, , , and °.
The drag time history indicates that, for °, °, and °, the maximum drag corresponds almost to the maximum lift and the minimum drag occurs near the minimum lift. At ° the trend of the variations is inverted; that is, increasing in the lift is related to a decrease in the drag and vice versa. The pitching moment oscillations follow the lift variations not the drag, because the lift variations amplitudes are bigger than the drag fluctuations. Thus, the lift is the determinative component of the pitching moment direction.
The computed mean values of the load coefficients are replotted in Figure 9 with corresponding range of the variations and the Strouhal numbers. Only the mean experimental load values are available for the modeled geometry. As previously observed, the calculated frequency of the fluctuations is similar for the vertical and horizontal forces. The amplitudes do, however, differ between the cases.
As mentioned, , is calculated based on the spoiler-iced airfoil projected height . The slight difference between the number of different angles of attack relates to different projected heights higher in °, though the frequency is almost the same. Then, there is a direct relationship between the air-projected height of the iced airfoil and the , while the amplitude of the oscillations does not correlate. Considering the separation bubble length as the length scale, the mentioned Strouhal number converts to ranging from 0.59 to 0.69 for different angles of attack. This corresponds to the shedding motion of the separation shear layer. It is consistent with the Strouhal number of regular mode for the same airfoil with a horn ice reported by Ansell [30]. Gurbacki [35] also reported similar range of (0.53–0.73) with a 3D horn ice shape.
Bigger “” value can be the source of more turbulent flow, that is, higher magnitude, at ° (Figure 6). It led to a larger trailing edge vortex shedding that was discussed before. It can be concluded that a higher turbulence intensity does not result in a higher amplitude of the oscillations.
(a) °
(b) °
(a) Pressure variation
(b) Schematic of the airfoil
(a)
(b)
(a) One instant of negative lift
(b) One instant of positive lift
The maximum amplitude of the lift fluctuations corresponds to zero angle of attack (Figure 9(a)). As mentioned, the lift force fluctuates between the negative and positive values at °. Flow behavior is analyzed in the following paragraphs for each case of negative, positive, and zero lift values. For the simulation of the flow at zero angle of attack, the vortices shapes at two time steps are shown in Figure 10. The time steps are selected with a half-period interval, in which the lift coefficient switches from a maximum positive to a minimum negative value.
Near the leading edge, a stagnation region arises, at which the flow is split into two bunches going to the upper and lower surfaces of the airfoil. They reach each other again near the trailing edge. Since the airfoil surface is covered thoroughly by the vortices, the streamlines on the upper and lower surfaces should traverse over the borders of the rotating flow regions. When there are two rotating bubbles on one side of the airfoil, the dead flow occupies a larger area. Consequently, the first active streamline should pass through a longer path. This path is considered from the splitting point upstream to the stagnation point downstream of the airfoil. Longer path leads to a higher velocity for the flow.
When the single bubble is on the upper surface (Figure 10(a)), the flow is slower at that side. This results in a higher pressure compared to the lower surface. The consequent higher pressure (less suction) on the upper surface justifies the negative lift. The reverse process happens in the case of positive lift (Figure 10(b)). The bubble configuration shifts and the stagnation point moves resulting in a shorter path for the lower streamline at another instant.
With similar arguments, the lift force becomes positive. These variations repeat periodically inducing a periodic variation of the lift force (Figure 8).
Previously, in this paper, the flow path was described regarding the pressure distribution and the approximate form of the separated region. Figure 10 indicates that the prescribed oval shape of the separated region may consist of more than one vortex.
Moreover, there are some fluctuations in the flow velocity downstream of the airfoil due to the shape and the location of the upstream vortices; that is, the vortices movement at the separation region affects the downstream flow where they are shedding. A similar period is obtained far downstream; that is, the fluctuations in the pressure and the velocity propagate at least four chord lengths downstream the iced airfoil which is modeled in this simulation. So, in practice, it can cause a periodic inlet velocity upstream of a neighbor wind turbine at wind farms. Therefore, for a wind farm in which the distance between the wind turbines is designed to be 3–10 rotor diameters, the icing effect should be considered as it may make time-varying instabilities downstream of the rotor.
Figure 11 shows streamlines in the proximity of the spoiler-ice, highlighting a multitude of vortices. It is an instant at which almost no lift is generated. The flow upstream of the spoilers forms a bubble in the corner between the spoiler-ice and the airfoil surface.
With the application of a RANS turbulence model these vortices are captured at the upstream of both upper and lower spoilers. The flow divides at the upstream face of the spoiler-ice and part of the flow is diverted to the upper spoiler while the remaining flow goes over the tip of the lower spoiler-ice and forms shear layers. These shear layers play a significant role in the generation of the downstream primary and secondary vortices as these vortices are clearly seen.
As the low speed rotating flows fill the region downstream of the spoilers, the airfoil surface is covered thoroughly by the vortices. Thus, the axial flow which was supposed to act on the airfoil surface to generate the lift force does not meet the airfoil surface at all. With the existing ice horn height, the designed curvature of the NACA0011 does not create any lift, as the flow is following the path formed at the boundary of the vortices instead.
Furthermore, although the integrated value of the vertical force on the iced airfoil is zero at the illustrated moment, there is a shedding vortex that initiates upward from the trailing edge (not shown here). Therefore, it is observed that even with a zero lift there is no symmetrical load distribution on the airfoil, and the downstream flow is oscillating.
As described in the literature, shear layer flapping is typically observed at very low frequencies, when compared to the characteristic frequencies of most other flow phenomena. A lower frequency cycle beside the main oscillations was observed when considering a long-time history. The related frequency is about 1/7 of the regular mode frequency, with the value giving at ° that is in consistent with the range of values reported in the literature for the simulated glaze ices [21, 26, 30].
4. Conclusions
The main objective of this paper is to estimate the dynamic loads acting over an iced airfoil, as well as studying the structure and dynamics of the turbulent separation bubbles. Ice profile is simulated with the help of spoilers. The simulated separated flow over the sharp spoilers can be considered as a worst test case of load loss due to the icing. It is shown that a glaze ice effect is not limited to a decrease in the lift; it also imposes some dynamic forces that should be considered in a wind farm, for instance.
Aerodynamic performance coefficients and pressure profiles are calculated and compared with the available measurements for a chord Reynolds number of 1.83 × 106 [9, 15]. The separation bubbles formed on an airfoil behind a simulated leading edge glaze ice accretion are studied by URANS. The details of the flow field and the vortex shapes at different angles of attack are investigated based on the mean pressure distributions as well as the instantaneous streamlines. The main conclusions can be summarized as follows:(i)The numerical results confirm both time-dependent phenomena observed in previous similar measurements [30]: a low-frequency mode, with a Strouhal number –0.02, and a higher frequency mode with a Strouhal number –0.69. The higher frequency motion has the same characteristics as the shedding mode and the lower one corresponds to the flapping mode.(ii)A greater pitching moment is a consequence of glaze ice beside the decrease of lift. Therefore, a higher angle of attack does not seem to be a good choice of operation in the presence of horn ice as changing the angle of attack to compensate lift loss can increase the pitching moment at the same time. Furthermore, the drag increases considerably.(iii)It is seen that the frequencies of the oscillations are almost the same for all angles of attack. There is a direct relationship between the iced airfoil projected height and the Strouhal number , while the amplitude of the oscillations does not correlate with that.(iv)A higher angle of attack leads to a higher turbulence intensity in the flow field, as the airfoil projected height increases.(v)The downstream oscillations propagate further downstream at higher angle of attacks.(vi)The drag variations in time are much smaller than the lift variation. It means that the vortices distribution is shifting more frequently between the upper and lower surfaces affecting mostly the vertical loads. However, both drag and lift oscillate at a similar frequency.(vii)Considering the load time histories of the simulated cases, maximum lift and maximum drag occur almost at the same time and the variations are in the same direction, while they are reverse at the negative angle of attack.(viii)The separated bubbles on both upper and lower surfaces are appearing closer to the leading edge compared to the experiments. In the experiments, the oval shape bubbles seem to be thinner near the extremities, while URANS modeling results in the vortices extending wider over the surfaces with small variation in the suction magnitude.(ix)A vortex arises close to the trailing edge at 10% of the chord length. The related vortex shedding is the source of the fluctuations downstream. Reaching more suction in that vortex leads to more fluctuations in the flow downstream.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this paper.