J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
Contents lists available at SciVerse ScienceDirect
Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia
Scaled ice accretion experiments on a rotating wind turbine blade Yiqiang Han n, Jose Palacios, Sven Schmitz The Pennsylvania State University, Department of Aerospace Engineering, University Park, PA 16802, United States
a r t i c l e i n f o
abstract
Article history: Received 27 October 2011 Received in revised form 18 May 2012 Accepted 1 June 2012
This work presents ice accretion tests of a model wind turbine blade mounted on a rotor test stand. The model test blade was scaled to reproduce local flow and icing conditions at the 95% radial station of the NREL Phase VI Rotor. A modified Ruff scaling method was implemented to scale atmospheric icing conditions that could be generated in a laboratory facility. Scaling laws allowed for reduced testing time as well as a reduced chord length of the blades compared to on-site icing events were implemented. Classical Blade Element Momentum Theory (BEMT) was utilized to correlate inflow conditions at the test facility to actual operating conditions of the NREL Phase VI Rotor. Experimentally obtained rime ice shapes are compared to LEWICE predictions and are used to validate the capability of the facility to reproduce representative icing conditions. The rime ice results matched with predictions to within 1% of both ice thickness and ice extent along the blade surfaces. The effects of angle-of-attack, temperature, liquid water content, and icing time on the final ice shapes were investigated in a parametric study. Experimentally obtained glaze ice shapes are used to demonstrate the severity of the icing events on wind turbine torque. The present work comprises testing procedures and experimental data that can be used for future efforts in ice accretion modeling and testing. & 2012 Elsevier Ltd. All rights reserved.
Keywords: Wind turbine icing Rotating test Scaling method Performance degradation
1. Introduction 1.1. Power loss due to wind turbine blade icing A large portion of the global wind resource exists in cold climate regions, where wind turbine performance is largely dependent on icing severity. Wind turbine icing results in blade aerodynamic performance degradation, increased blade fatigue, and safety concerns due to ice shedding along with rotor imbalance. Safety concerns often lead to turbine shut down during parts of the winter season, thus affecting the Annual Energy Production (AEP) (Homola et al., 2010). For example, the AEP loss due to icing at the ¨ Alpine Test Site Gutsch (2350 m above sea level) was estimated to be 23% based on field measurements (Barber et al., 2011). At Nordland in Norway (770 m above sea level), the Weather Research and Forecast (WRF) code was used to study wind power loss due to icing (Byrkjedal and Vindteknikk, 2009). The estimated icing time was 1000 hours per year, and the corresponding AEP loss ranged between 14% and 28%. A computer simulation of an icing event and the resulting power loss has also been conducted for the NREL 5 MW baseline wind turbine (Homola et al., 2012). The estimated power loss due to icing was found to be between 24% and 27%.
n
Corresponding author. Tel.: þ1 81 4867 2982. E-mail address:
[email protected] (Y. Han).
0167-6105/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jweia.2012.06.001
Although the possible loss of power production due to icing has been well recognized, wind energy in cold climate regions still remains desirable. These regions have in general high winds and the air density could be 27% higher at 30 1C compared to that at þ 35 1C. Mountainous wind turbine sites in the centraland north-eastern U.S., northern Europe, and some parts of Asia are rich in wind resource; whereas, in contrast, those sites are subject to adverse icing environments that may last for several months. During an atmospheric icing event, when liquid water still exists in ground-level clouds at temperatures below freezing, rotating wind turbine blades start to accrete ice on the leading edge. The associated surface alteration of sectional airfoils leads to instant performance degradation, which is reflected in reduced lift-to-drag ratios. Airfoil performance degradation due to icing was quantified in a wind tunnel experiment (Jasinski et al., 1998) at the University of Illinois at Urbana-Champaign (UIUC) on a S809 airfoil (Somers, 1997) equipped with fabricated ice shapes obtained by NASA LEWICE simulations. For a typical wind turbine rime ice shape and an ice thickness of only 2.5% of the blade chord, the resultant loss in power coefficient, CP, was identified to be 14.5% under a variable operational speed mode using a computer model based on Blade Element Momentum Theory (BEMT). At the primary monitoring location (75% r/R), a 168.2% increase in drag coefficient, Cd, was associated with the significant loss in power coefficient, CP. In light of the significant performance degradation caused by icing events, ice accretion experiments on wind turbine airfoils are
56
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
desirable. A wealth of literature on airfoil performance degradation under aircraft icing conditions is available (Flemming and Lednicer, 1986; Reinert and Aubert, 2011; Han et al., 2011). An equivalent database is required but to date does not exist for wind turbines under representative atmospheric ground icing conditions. Systematic performance evaluation of wind turbine airfoils under icing conditions needs to be conducted to provide a valuable database for modeling of wind turbine operation in cold climate regions. The major obstacle of the on-site icing testing is that data acquisition of ice accretion events at wind turbine sites is a difficult task due to uncontrolled icing conditions. Obtaining ice shapes directly on wind turbine blades is only possible through photographs and is hence of limited accuracy. Controlled environments such as icing wind tunnels have provided ice shapes on generic airfoils which have been correlated to 2D ice prediction codes (Shin et al., 1994). A very limited number of experiments have been performed on a rotating test stand (for example, Flemming and Lednicer, 1986). Test data on torque increase are not available from stationary icing wind tunnel tests. An ideal 2D test section in an icing wind tunnel is desirable to compare ice accretion for fixed wing aircraft applications. For wind turbine applications, on the contrary, rotational effects on ice accretion need to be further understood. A rotating test stand provides an opportunity to study these effects. Another concern related to ice accretion tests is that, in a natural icing event, the icing duration could be hours or even several days. It is impractical to test the same sized wind turbine airfoil under the same icing conditions for such a long time in an icing wind tunnel or other research facility. An icing scaling method can be very beneficial in accelerating ice accretion tests by using a smaller airfoil chord, higher free stream velocity, and by providing the incoming flow with higher water content. Such a scaling method needs to go through rigorous examination to ensure that the same normalized ice shapes can be obtained. 1.2. Ice shape categorization Depending on the icing conditions, ice shapes accreting on rotating blades can be categorized into three groups: glaze ice, rime ice, and mixed ice (Heinrich et al., 1991). Glaze ice is formed by water droplets that impact and freeze as they travel along the blade surface. It occurs at relatively warm temperatures (0–5 1C) with Liquid Water Contents (LWC) larger than 0.2 g/m3 and at high air speeds. Rime ice, on the other hand, is formed when water droplets freeze immediately upon impingement. This type of ice occurs mainly at colder temperatures ( 10 1C or lower) with low LWC (below 0.2 g/m3) (Belte, 1987). Mixed ice is categorized as a blend between glaze- and rime ice, and as such it exhibits both the round leading edge rime ice shape and the protruding ice feathers at the back of the ice shape that resemble glaze ice. Sample ice shapes for all three ice types are shown in
Fig. 1 on a generic NACA 0012 airfoil from results of previous icing tests (Han et al., 2011). In general, glaze ice introduces a more severe drag penalty than rime ice, since it creates an irregular airfoil surface and often has two protruded horns that considerably decrease the lift-todrag ratio. During rime ice accretion, water droplets freeze upon impact, thus eliminating complicated 3D interference of surface running water from other locations along the blade span. Rime ice is often chosen for model validation due to its simpler accretion mechanism. 1.3. Objectives The objectives of this research are: (1) to generate a practical icing condition matrix applying existent ice scaling methods, (2) to obtain ice shapes on a rotating wind turbine blade, (3) to validate the capability of the facility to reproduce wind turbine representative rime ice shapes, and (4) to conduct a parametric study on icing parameters and ice severity for both rime and glaze ice shapes.
2. Test apparatus and methods 2.1. Icing test facility The Adverse Environment Rotor Test Stand (AERTS) was designed and constructed at the V/STOL Research Center of Excellence at the Pennsylvania State University (Palacios et al., 2012). The purpose of this facility is to advance the understanding of the fundamental physics of ice accretion and to evaluate ice protection systems. Fig. 2 shows the AERTS facility with a set of 1.372 m (4.5 ft) radius S809 wind turbine blades. The AERTS test section is formed by a cold chamber with dimensions of 6 m 6 m 3.5 m. A 24 signal channel/24 power channel slip ring is used for signal transmission to the control console and to power ice protection systems located on the rotating blades. During an ice accretion test, the rotor is driven by an 89.5 kW (120 Hp) motor to maintain a specified rotational speed up to 1000 RPM. The icing cloud is generated by standard icing spray nozzles that were donated by the NASA Glenn Icing Research Tunnel (IRT). When the blade encounters the sprayed water droplets at low ambient temperatures, ice starts to accrete at the blade leading edge. During icing tests, forces and moments acting on the test stand are measured by a 6-axis load cell built in the bell-housing. The shaft torque required to maintain the specified rotational speed is measured by a sensor placed in-line with the drive shaft. After the desired icing conditions are maintained for the duration of the icing event, resultant ice shapes can be measured at different monitor locations along the blade span. Measured ice shapes can serve as validation cases for
Fig. 1. Three types of accreted ice on NACA 0012 airfoil from previous AERTS tests.
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
57
measuring the ice thickness at the stagnation line of the test blade. A LWC determination code utilizing a set of ice scaling equations was used to calculate the experimental LWC based on stagnation ice thickness and other icing parameters, such as velocity, temperature, MVD, etc. (Palacios et al., 2012).
2.2. Scaling icing conditions—modified Ruff method
Fig. 2. AERTS test stand with S809 test blade.
ice prediction models and other CFD studies; measured force and torque data can be used to evaluate blade performance degradation. In addition to ice accretion tests, the facility allows shedding tests used for the evaluation of blade coatings. Anti-icing/de-icing system performance tests can also be conducted. The icing conditions that can be controlled and measured in the facility include: temperature, icing time, rotor speed, Liquid Water Content (LWC) and Median Volume Diameter (MVD). The temperature can be controlled in the range 0 to 25 1C and kept constant during a typical 30 min ice accretion test. To match the operational speed of a wind turbine, a rotor RPM ranging between 300 and 500 must be maintained. These RPMs correspond to a wind turbine tip speed varying between 43 and 72 m/s. To physically characterize an icing cloud, two icing parameters, known as Liquid Water Content (LWC) and Median Volume Diameter (MVD) are used. The LWC is the characteristic waterto-air concentration in a two phase flow (liquid and gas). The MVD is a characteristic number that describes the average water droplet size in an icing cloud. The value of MVD is obtained by assuming that half of the liquid water volume in the icing cloud is contained in drops larger and half in drops smaller than the quoted value. It has been demonstrated that the drop motion in the flow field can be modeled using a single value of MVD instead of a certain statistical size distribution (Finstad et al., 1988). There are very limited recorded LWC and MVD datasets available that were directly measured at wind turbine sites. At Murdochville, Canada (900 m above sea level), LWC and MVD in two observed icing events have been reported (Hochart et al., 2008). Recorded LWCs were 0.218 g/m3 and 0.242 g/m3 with MVDs of 38.3 mm and 40.5 mm, respectively. Another systematically measured set of ¨ (719 m above sea icing conditions was recorded at Mount Yllas level) (Nygaard et al., 2011). The measured LWC ranged from 0.09 to 0.43 g/m3 and MVD between 12.1 and 19.9 mm. Based on this information, the icing envelope for wind turbine ice testing is determined to be: (i) LWC 0.1 g/m3 and MVD o20 mm for rime ice cases; and (ii) 0.2 g/m3 oLWC o1 g/m3 and MVDo40 mm for glaze ice cases. The direct measurement of MVD and LWC using Phase Doppler Particle Analyzer (PDPA) and Forward Scattering Spectrometer Probe (FSSP) at the rotor plane is not appropriate in AERTS, since the high centrifugal force inherent in the rotating frame may pose potential harm to the equipment. Instead, the MVD is directly controlled by adjusting the input air–water pressure ratio to the spraying nozzles. The MVD can be readily calculated from the NASA calibrated nozzles. The LWC was measured for every test by
Due to size limitations in icing facilities, an icing scaling method needs to be used with the aim of providing identical non-dimensional ice shapes to both the sub-scale experimental results and the actual wind turbine. The icing scaling method is needed, since a given icing facility may not be able to achieve a certain range of icing conditions in terms of velocity, temperature, model geometry, and icing cloud (Makkonen and Oleskiw, 1997). By applying a scaling law to a sub-scale model, the required icing time can also be reduced. An accelerated icing experiment can be conducted instead of testing full-scale models for long periods of time. This accelerated testing is very favorable for wind turbine icing research, since the typical icing duration in nature is quite long, usually hours or even days. In general, models with smaller chord lengths are tested in icing facilities. However, an identical normalized ice shape can still be obtained by matching the flowfield and icing parameters at a particular spanwise location along the wind turbine blade. Local flow conditions for the actual operating blade are calculated by BEMT. The local inflow conditions consisting of air speed and angle of attack (AOA) are the same at the desired location of the actual rotor and the test blade in the AERTS facility. A modified Ruff method (Ruff, 1986) is used to scale the icing conditions. The modified Ruff method is a validated NASA scaling method that was originally designed for aircraft icing tests in wind tunnels (Anderson, 2004). By matching the characteristic numbers of six similarities between two different models, an ice scaling similarity can be established. A flow chart of the scaling method is shown in Fig. 3. The normalized stagnation line ice thickness at the blade leading edge can be calculated during the first three similarity analyses in Fig. 3. The same non-dimensional iced profile can be achieved by matching
D d
¼ n0 Ac b0
ð1Þ
where D is the stagnation line ice thickness and d is the characteristic length, which is twice the blade leading edge radius in this case. The right-hand-side of Eq. (1) is a product of three important icing parameters obtained from the first three analyses in the scaling method. Here the collection efficiency, b0, with respect to the droplet trajectory similarity, was defined to illustrate the fraction of incoming water content that actually impacts the monitoring control volume (Langmuir and Blodgett, 1946). The expression for b0 is
b0 ¼
1:40ðK 0 ð1=8ÞÞ0:84 1þ 1:40ðK 0 ð1=8ÞÞ0:84
ð2Þ
where the subscript, 0, denotes that it is calculated at the stagnation line. It is assumed that there is no incoming interfering water into the control volume at this location. K0 is a modified inertia parameter. It is a function of MVD, impact velocity, air viscosity, air density, and water density defined by Langmuir and Blodgett. During the water catch similarity analysis, the accumulation parameter, Ac, was defined to show normalized maximum local
58
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
Fig. 3. Flow chart of scaling method.
ice thickness Ac ¼
LWC V t ri d
ð3Þ
where V is the free stream velocity, t is the icing time, ri is the density of ice, and d is the characteristic model dimension, which is again twice the leading edge radius for symmetric airfoils. The freezing fraction, n0, is introduced next to denote the ratio of impinging water that freezes within a control volume. This term was first introduced by Messinger (1953) and later further developed by Ruff (1986). It can be expressed as follows: C p,ws y n0 ¼ ð4Þ fþ Lf b where several characteristic energy coefficients are comprised: Cp,ws is the specific heat of water on the model surface; Lf is the latent heat of freezing water; f is the drop energy transfer parameter; y is the air energy transfer parameter, and b is the relative heat factor. 2.3. Geometric specifications of wind turbine test blade The AERTS facility, shown in Fig. 2, is a rotor test stand operating in a controllable icing cloud. It was found in previous research that conditions at the 99% radial station of an AERTS blade are twodimensional to a very good approximation (Han et al., 2011). The ground effect and wall effect are considered secondary effects for the rotor blade inflow for lightly loaded rotors. Due to ease of manufacture and previous experience with inflow distribution inside the facility, a rectangular and untwisted wind turbine blade equipped with the S809 airfoil was designed and fabricated to be tested in the AERTS facility. The blade length is 1.016 m (40 in.), the blade chord is 0.267 m (10.5 in.), and the aerodynamic part of the blade is attached to a root spar section such that the actual blade radius is 1.372 m (54 in.); see Fig. 4. Ice shapes and thickness were monitored at 6 locations along the blade span (51.9%, 61.5% 71.2%, 80.8%, 90.4% and 99%). The primary monitoring location was at the 99% radius. Scaled icing conditions were selected to match at this location, and ice shapes were recorded at that location for comparison with wind tunnel data and ice accretion simulations. The reference case for all ice accretion tests of this work is the NREL Phase VI Rotor (Hand et al., 2001). The NREL Phase VI Rotor is a two-bladed and stall-controlled wind turbine, whose twisted
Fig. 4. AERTS S809 test blade.
Table 1 Basic specifications of the NREL Phase VI rotor. NREL Phase VI rotor Radial location (r/R) Local radius (m) Chord (m) Tip pitch angle (deg.) Rotational speed (RPM)
95% 4.778 0.3762 31 72
and tapered blades are exclusively equipped with the S809 airfoil. Some rotor specifications are comprised in Table 1. The monitoring location on the target rotor blade is at the 95% radial location. The reason for choosing this location is that the ice load at the tip is of highest interest, as the bulk of produced rotor torque and hence loss of rotor torque due to icing occurs in this region. In another ice modeling research using the same NREL VI turbine rotor (Fu and Farzaneh, 2010), it was concluded that the ice load and associated loss in rotor power shows an increase towards the blade tip, while little ice can be observed at the blade root. Another research at a wind turbine test site has reported that icing in the blade tip region, 95–100% blade span, has the most pronounced effect on the wind turbine performance (Barber et al., 2011). In addition, the following characteristics at the 95% radial
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
station are favorable to the icing scaling method, i.e. the local velocity of the ice encountered at the tip (36–48 m/s, corresponding RPM 72–96) is representative of a typical ice accretion speed at an icing test facility. The chord at the blade tip is also smaller than that at the inboard area, which is 0.3762 m. The chord of the test blade in the AERTS facility is 0.2667 m, which gives a scaling ratio of 1.4. A scaling ratio less than 2 is favorable for the reproduction of practical ice shapes. To provide aerodynamic inflow conditions for the scaling calculation, Blade Element Momentum Theory (BEMT) was used to find the required sectional air speed (V) and AOA (a) at the monitoring locations in addition to the corresponding collective pitch of the test blade. In classical BEMT, the sectional air speed, V, is related to the wind speed, V0, via V¼
V 0 ð1aÞ Orð1 þa0 Þ ¼ sin F cos F
ð6Þ 0
where a is the axial induction factor, a is the angular induction factor, and
F ¼ a þ Db95%-Tip þ bTip
ð7Þ
where F is the local blade flow angle and bTip is the blade tip pitch angle. For the NREL Phase VI Rotor, the offset in local blade pitch between the 95% radial station and the blade tip, Db95%-Tip in Eq. (7), is 0.3471 The induction factors, a and a0 , are functions of the local blade flow angle, F, and the local AOA, a. They can be described as (Manwell et al., 2010) ( )1 4FðFÞsin2 ðFÞ a ¼ 1þ 0 ð8Þ s ½C l cos ðFÞ þ C d sin ðFÞ a0 ¼ 1 þ
1 4FðFÞsinðFÞcosðFÞ s0 ½C l sinðFÞC d cosðFÞ
ð9Þ
Table 2 NREL Phase VI rotor—operating conditions considered for AERTS tests. Operating Case #
V0 RPM bTip (m/s) (deg.)
Remarks
V (m/s) AOA (deg.)
1 2
7.5 7.5
72 72
1.0 10.75
Actual NREL Phase VI Rotor
36.1 36.1
6 0
3 4 5 6
6.5 6.0 9.5 4.0
96 84 84 84
6.5 0.0 0.0 0.0
Variable-Speed NREL Phase VI Rotor
47.9 42.0 42.0 42.0
0 4 8 2
59
In Eqs. (8) and (9), F(F) is the total loss factor that is, in general, the product of Prandtl’s tip-loss factor and an adjusted root loss factor (Glauert, 1935; Manwell et al., 2010). s0 is the local rotor solidity. The sectional airfoil lift- and drag-coefficients, Cl and Cd, are functions of the local blade flow angle, F, and AOA, a. Given the air speed, V, and the rotor speed, O, Eqs. (6)–(9) can be used to find the unknowns, i.e. the wind speed, V0, the blade tip pitch angle, bTip, and the local AOA, a. As Eqs. (6)–(9) describe a non-linear system of equations, the experienced blade designer can find different combinations of V0 and bTip for a desired a; likewise a desired sectional air speed, V, can be obtained by various combinations of V0 and bTip with a being a result of the analysis. In this work, a total of 6 operating conditions were considered for the NREL Phase VI Rotor. They are summarized in Table 2. An in-house wind turbine blade design and analysis code, XTurb-PSU (Schmitz, 2011), was used to solve Eqs. (6)–(9). The XTurb-PSU code can be run in two modes, i.e. (i) BEMT mode rooted in NREL’s AeroDyn model (Moriarty and Hansen, 2005), and (ii) Lifting-Line mode based on a prescribed vortex method (Chattot, 2005). The BEMT mode was used for all analyses presented in this work. Fig. 5 illustrates computed AOAs, a, and local air speed, V, along the blade span of the NREL Phase VI Rotor. Values at the 95% radial station (highlighted with rectangles) correspond to the inflow conditions tested in the AERTS facility. In the BEMT analysis, the wind speed, V0, and blade tip pitch angle, bTip, were adjusted manually such that the air speed, V, and the local AOA, a, were within 2% of the desired values.
2.4. Ice accretion test matrix In total, 17 ice accretion tests were conducted during this research, see Table 3. For all of the experimental cloud, LWCs were less than 1.5 g/m3, which is a representative range encountered by a wind turbine during the winter season. For comparison of experimental rime ice shapes with predictions, LWCs were maintained very low (0.08 g/m3 and 0.05 g/m3) in order to ensure rime ice accretion with an icing time of 30 min for a single test. After applying the scaling laws, the icing time required in AERTS to reproduce a wind turbine ice shape is approximately half the time of the icing duration encountered by the wind turbine. The maximum icing time tested during this research was 162 min (AERTS case number 20–22). The tested icing times are representative of an average intermittent rime icing event experienced by wind turbines in the field.
Fig. 5. NREL Phase VI Rotor—BEMT analysis and operating cases at 95% r/R.
60
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
Table 3 Test matrix—AERTS S809 blade. Test Conditions for AERTS S809 Rotor
Scaling Conditions for NREL Phase VI Rotor
AERTS case #
MVD (lm)
LWC (g/m3)
Temp. (1C)
Vel. (m/s)
RPM
Time (min)
AOA (deg.)
MVD (lm)
LWC (g/m3)
Rime 1 2 3 4 5 6 7 8 9
ice 4 15 16 17 18 19 20 21 22
20 20 20 20 20 20 20 20 20
0.08 0.08 0.08 0.08 0.05 0.05 0.05 0.05 0.05
7 9 4.5 7 10 9 9 9 9
50 50 50 50 50 50 50 50 50
350 350 350 350 350 350 350 350 350
30 30 30 30 30 30 30 60 90
2 4 2 4 4 8 4 4 4
26.6 26.6 26.6 26.6 26.6 26.6 26.6 26.6 26.6
0.076 0.076 0.077 0.076 0.047 0.047 0.047 0.047 0.047
Glaze 10 11 12 13 14 15 16 17
ice 1 2 9 12 14 23 24 25
20 30 20 20 20 20 20 20
1.2 1.2 0.2 0.3 0.4 0.3 0.3 0.3
16 8.3 16.4 11.4 8.5 8.5 9 4.5
57 57 43 43 43 43 43 43
400 400 300 300 300 300 300 300
6 4 4 4 25 4 12 12
0 0 0 0 0 6 6 6
26.6 39.9 26.6 26.6 26.6 26.6 26.6 26.6
1.12 1.13 0.19 0.28 0.37 0.28 0.28 0.28
The glaze ice shapes on the AERTS scaled rotor were obtained under relatively high LWCs (ranging from 0.2 to 1.2 g/m3) and for shorter icing times (4–25 min). The reduced time during glaze ice testing was used to prevent the potential of ice shedding events and the consequent rotor imbalance. The corresponding glaze icing times for the full-scale NREL Phase VI rotor thus range between 7 min and 45 min. Such conditions represent severe icing events that wind turbines may experience during intermittent adverse weather.
3. Ice accretion experimental results In order to validate the capability of the AERTS facility in reproducing wind turbine ice shapes, 9 rime ice shapes were compared to predictions obtained by the LEWICE 2D version 3.2 (Wright, 2008). The LEWICE ice prediction code was developed originally by the NASA Glenn Research Center (Ruff and Berkowitz, 1990) and has been validated against extensive datasets from experiments conducted at the Icing Research Tunnel (IRT) at NASA Glenn (Wright, 2005). It has been prevalently used for decades by the aircraft industry for safety certification. It has also been introduced to the wind industry for ice prediction and performance evaluation on wind turbines (Jasinski et al., 1998; Barber et al., 2011). Other ice prediction codes have been developed, such as TURBICE (Makkonen et al., 2001) and FENSAP-ICE (Bourgault et al., 2000), however this work used the LEWICE 2D prediction code exclusively due to availability and prior experience with rotorcraft validation studies in AERTS. The LEWICE code is known to be capable of accurately predicting dry ice growth (rime ice). The wet growth (glaze ice), on the other hand, is difficult to predict since inherent complexity results from surface running water after water droplet impact. It has been observed by researchers that LEWICE cannot accurately predict ice shapes for certain glaze ice conditions (Wright, 2008). This is mainly caused by the complicated surface energy exchange under the wet growth condition. The difference between glaze ice prediction and NASA icing experiments has also been observed by researchers using other ice prediction codes, such as TURBICE
Temp. (1C)
6.9 8.9 4.4 6.9 9.9 8.9 8.9 8.9 8.9
15.9 8.2 16.3 11.3 8.4 8.4 8.9 4.4
Vel. (m/s)
RPM
Time (min)
AOA (deg.)
Table 2 Operating case #
42 42 42 42 42 42 42 42 42
84 84 84 84 84 84 84 84 84
53.8 54.1 53 53.8 54.2 54.1 54.1 108.2 162.3
2 4 2 4 4 8 4 4 4
6 4 6 4 4 5 4 4 4
47.9 47.9 36.1 36.1 36.1 36.1 36.1 36.1
96 96 72 72 72 72 72 72
10.9 7.1 7.3 7.3 45.4 7.3 21.8 21.5
0 0 0 0 0 6 6 6
3 3 2 2 2 1 1 1
(Makkonen et al., 2001) and THERMICE (Tran et al., 1994). So far, the surface roughness and the corresponding heat transfer balance on the wet surface at the initial stage of icing is still challenging to predict accurately, and a limited number of experimental datasets exists for validation of available prediction codes. Also, the centrifugal force inherent in the rotating frame tends to drive the surface running water outboard (Han et al., 2011). This kind of 3D effect on the ice shape is a unique feature of the rotating ice accretion under glaze icing conditions, which cannot be captured by 2D ice prediction codes. The 8 tests on glaze ice shapes, see Table 3, were therefore not compared to LEWICE computations. On the contrary, these ice shapes are thought to serve as an experimental database to the wind energy community for subsequent efforts on glaze ice accretion modeling. For the glaze ice shapes, the measured rotor torque was analyzed during the ice severity study.
3.1. Experimental ice reproduction validation and parametric study of icing parameters Nine (9) rime ice shapes obtained in the AERTS facility are compared to LEWICE predictions. A variety of icing conditions were tested, see Table 3. The 6 ice shapes in Figs. 6 and 7 show the 30 minute test results for different LWC values (0.08 and 0.05 g/m3), temperatures (4.5 1C, 7 1C, and 9 1C), and AOAs (21, 41, and 81). Ice thicknesses at the leading edge stagnation point of these ice shapes amount to less than 2 percent of the blade chord, thus corresponding to light rime icing. The ice thicknesses and overall rime ice shapes match very well between the prediction code and the experiments. This indicates the capability of the AERTS facility in reproducing and controlling icing conditions around a rotor test stand. The effects of various icing parameters are also examined in this section. Ice shapes with single-value variations of AOA, temperature, LWC, and icing time are compared, while other conditions are maintained ceteris paribus. During ice accretion testing, the MVD was maintained at 20 mm for all tests. Test velocities varied to provide three conditions: 43 m/s (300 RPM), 50 m/s (350 RPM) and
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
61
Fig. 6. Rime ice shapes obtained at varying temperature and angle-of-attack—1: (a) AERTS case 4, (b) AERTS case 15 and (c) AERTS case 16.
57 m/s (400 RPM). In this research, the effects of MVD and velocity on the final ice shape were not studied.
correctly is also important for the accurate prediction of the blade’s performance degradation due to icing.
3.1.1. Effect of Angle-of-attack (AOA) The effect of AOA on the ice shapes can be identified in Fig. 8. Small AOAs provided no significant difference in accreted ice shapes (see Fig. 8(a)). Changes become more evident as the AOA is further increased, see Fig. 8(b). Since wind turbines typically operate at AOAs between 41 and 101, an accurate matching of the local AOA is critical to ice accretion modeling. In addition, the AOA will vary with wind conditions during an icing event. Capturing these AOA changes
3.1.2. Effect of temperature The growth of two main ice shape categories, rime- and glazeice, are dominated by icing temperature. For rime ice cases shown in Fig. 9, the effect of temperature is not significant. It is worth mentioning that these rime ice cases were conducted at low LWC and long icing time conditions, see Table 3. The little amount of incoming water droplets in the test chamber resulted in well supercooled water droplets that froze immediately upon impingement with the blade. The instantaneous water freezing in the rime regime
62
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
Fig. 7. Rime ice shapes obtained at varying temperature and angle-of-attack—2: (a) AERTS case 17, (b) AERTS case 18 and (c) AERTS case 19.
eliminates water run-back. Comparison of Fig. 9(a) and (b) reveals that as temperature became colder, discrepancies in accreted ice shapes became smaller. This is attributed to the fact that less water splashed or traveled in the chordwise and spanwise directions, which suggests that the rotational effect is not dominant in rime ice cases. In general, temperature has a very small effect on wind turbine rime ice shapes. For glaze ice shapes in Fig. 10, the situation is different. Since more water droplets (LWC¼0.3 g/m3) impact onto the airfoil at warmer temperature (T¼ 4.5 1C), surface running water effects are likely to occur. It can be seen from Fig. 10 that the leading edge ice thickness and upper surface ice impingement limit matched between the two test cases. This indicates that the total amount of incoming water is the same and independent of the temperature difference. On the lower surface, however, an apparent difference
can be seen in the surface impingement limits. This difference is attributed to surface running water that originated from the inboard part of the blade span and moved outboard under the action of centrifugal forces. The ice growth limit of the 4.5 1C case is further downstream than that of the 9 1C case. This spanwise surface running water condition at the blade tip is less significant for helicopters, since the relatively high rotational speeds of helicopter blades (200–400 RPM) can generate sufficient centrifugal force to shed off the ice at the blade tip. Furthermore, kinetic heating in that area due to the high local speeds (up to 750 ft/s) also prevent large amounts of ice accretion at that location. This is different for wind turbines, since most of the wind turbines operate at low rotational speed (8–20 RPM), and it is known that ice accreted at wind turbine blade tips generates the largest blade loads, ice shedding concerns, and contributes greatly to performance degradation (Fu and
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
Farzaneh, 2010). For wind turbines, the surface running water phenomenon adds to the turbine performance degradation and thus to overall loss in energy production. Modeling this type of ice accretion requires three-dimensional solvers with accurate water run-back and splashing models. Preliminary experimental data obtained in this work can serve as a future validation database for full three-dimensional ice accretion models. 3.1.3. Effect of LWC The effect of LWC is shown in Fig. 11. As expected, the ice thickness varies approximately linearly with changes in LWC. This linear increase is attributed to the relationship between the mass flow rate of incoming water droplet and the LWC, i.e. _ ¼ LWC V b0 m
ð10Þ
where V is the free-stream velocity and b0 is the collection efficiency (Langmuir and Blodgett, 1946). 3.1.4. Effect of time Two sets of ice shapes are presented in Fig. 12. As it can be seen from Fig. 12(a), the three tests under rime ice conditions show fairly good uniformity with icing time. The mass of frozen water droplets in the three ice layers (corresponding to the three icing times) are approximately the same. For glaze ice cases in Fig. 12(b), the 4 min test resulted in a smooth ice shape similar to the rime ice shapes. This kind of ice shape can still be accurately predicted by LEWICE despite being in the glaze regime. For the 12 min glaze test, the ice growth limit extended downstream towards the trailing edge. The ice surface became irregular as run-back effects became more severe, which resulted in higher aerodynamic penalties in terms of drag rise and associated rotor torque (ultimately affecting wind energy production). Prolonged glaze ice accretion on wind turbines can only be modeled using sophisticated three-dimensional ice accretion codes that include accurate run-back and splashing models. These data can assist in validating the many methods that are currently under development, such as performance degradation predictions using the 2D FENSAP-BEMT method (Homola et al., 2011), 3D Eulerian–Eulerian method with Fluent k–e solver (Fu and Farzaneh, 2010), 3D LES–Lagrangian method (Szasz and Fuchs, 2011), etc. 3.2. Performance degradation study Ice accretion leads to aerodynamic performance degradation in terms of decreased airfoil lift coefficient, Cl, and increased drag coefficient, Cd. The thrust produced by the AERTS test rotor and its corresponding change due to icing were fairly low; the typical thrust generated by the S809 rotor was 133.4 N (30 lbf) at 81 collective pitch angle. The drag rise results in increased torque requirements
63
to sustain a desired rotor RPM which is an indicator of the icing severity. This study of the performance degradation focused on increased rotor torque due to an increased blade drag coefficient. The torque at the driving shaft was monitored by a torque sensor and recorded by the data acquisition system in the AERTS facility. Besides the torque rise due to blade icing, other contributors to the measured rotor torque are friction in the transmission and bearings. Initial experiments to evaluate these torque losses were conducted prior to icing tests by measuring the required torque for the rotating hub without blades at rotational speeds between 300 RPM and 400 RPM. Transmission losses were quantified to be 16 N m. A BEMT rotor performance code was then used to estimate the required driving torque for the clean blade. The predicted torque data are the summation of the torque requirement for the clean blade and the torque loss in the rotor stand transmission. The torque requirements at the beginning and end of each test are summarized in Table 4. The corresponding icing conditions are also listed. Torque requirements for both rime- and glaze-ice cases were recorded during the tests. Ice accretion tests were categorized into three groups: (i) light-, (ii) moderate-, and (iii) severe-ice coverage. A comparison between torque requirements for all three categories is shown in Fig. 13. For light ice coverage, the performance losses were less than 5%. Moderate ice coverage resulted in an approximate 25% increase in required torque. Under severe ice coverage, the rotor required more than 70% torque to maintain the desired RPM. The total rotor torque data presented herein do not allow the extraction of local airfoil drag coefficients because of spanwise recirculation phenomena in the test chamber that are currently not being modeled and understood. The torque measurements illustrate the potential impact of ice accretion on performance degradation of a wind turbine operating in such adverse conditions. 3.2.1. Light ice coverage AERTS cases 15 and 17 are representative of light ice coverage (see Fig. 9(b)). Both cases were conducted under the same rime ice conditions, however at slightly different temperatures, see Table 3. It was shown though in Fig. 9(b) that temperature in this range has only a very small effect on the resultant ice shape. A comparison of the time-history of the rotor torque is shown in Fig. 14(a); the associated torque increments are shown in Fig. 13. The torque recorded during case 23, shown in Fig. 14(b), depicts a decrease in torque requirements to maintain a desired RPM. This effective increase in rotor performance due to ice accretion is attributed to the ice shape in Fig. 12(b) acting as a leading edge flap. The collective blade pitch angle for this case was 14.61 which resulted in local AOAs larger than 61 along the whole blade span, see Table 3. The cambered shape delays stall at the local airfoil section, thus increasing the rotor performance. This phenomenon was first recorded during icing tests at the Lewis Flight Propulsion
Table 4 Effect of icing on rotor torque—AERTS S809 blade. Vel. (m/s)
RPM
Time (min)
AOA (deg.)
Collective pitch (deg.)
Predicted baseline torque (N m)
Initial torque (N m)
Final torque (N m)
Torque increment (%)
Remarks
9 7 8.5
50 50 43
350 350 300
30 30 4
4 4 6
11 11 14.6
49.9 49.9 55.4
51.9 50.0 55.1
53.8 52.4 52.3
3.7% 4.8% 5.1%
Light Ice Coverage
0.2 0.3
16.4 11.4
43 43
300 300
4 4
0 0
1 1
19.6 19.6
20.3 21.0
25.2 26.9
24.1% 28.1%
0.4 1.2 1.2
8.5 16 8.3
43 57 57
300 400 400
25 6 4
0 0 0
1 1 1
19.6 20.4 20.4
19.2 17.4 19.5
32.8 37.9 38.2
70.8% 117.8% 95.9%
AERTS Case #
MVD (lm)
LWC (g/m3)
1 2 3
15 17 23
20 20 20
0.08 0.08 0.3
4 5
9 12
20 20
6 7 8
14 1 2
20 20 30
Temp. (1C)
Moderate Ice Coverage Severe Ice Coverage
64
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
Fig. 8. Effect of AOA on rime ice shapes: (a) LWC¼0.08 g/m3 and (b) LWC¼ 0.05 g/m3.
Fig. 9. Effect of temperature variation—rime ice (a) LWC¼ 0.08 g/m3 and (b) LWC¼ 0.08 g/m3.
Laboratory, now the NASA Glenn Research Center (Gray and von Glahn, 1958). It was claimed by the authors that rime ice formations at higher AOA often resemble depressed nose flaps; such ‘‘flaps’’ apparently assist airflow over the nose and reduce the extent of flow separation. This kind of leading edge flap phenomenon has also been observed and reported by other authors on a wind turbine airfoil (Jasinski et al., 1998). 3.2.2. Moderate ice coverage AERTS cases 9 and 12 are considered to be moderate rime ice shapes. These two cases are representative of extreme adverse weather conditions at low temperatures ( 11.4 1C and 16.4 1C) and high LWC (0.2 g/m3 and 0.3 g/m3). The increments of torque for these two cases are 24.1% and 28.1%, respectively, see also Fig. 15. Note that these tests only lasted for 4 min, which is equivalent to 7.3 min for a NREL Phase VI Rotor operating in the atmosphere. Such an icing condition could last for a prolonged time, and it is highly probable that the turbine would have to be shut down under these conditions.
Fig. 10. Effect of temperature variation–glaze ice (LWC ¼0.3 g/m3).
3.2.3. Severe ice coverage Three severe ice cases (AERTS cases 1, 2, and 14) are shown in Fig. 16. These are glaze ice shapes at high LWC (0.4 g/m3 and
1.2 g/m3). In general, these ice shapes are encountered during intermittent icing environments, where icing conditions can fluctuate vigorously (Heinrich et al., 1991). The torque measured
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
65
Fig. 13. Relative changes in rotor torque due to icing severity—AERTS S809 Blade.
Fig. 11. Effect of LWC variation—rime ice (LWC ¼ 0.05 g/m3 and 0.08 g/m3).
Fig. 14. Torque measurement (Light Ice Coverage)—AERTS S809 blade: (a) rime ice (LWC ¼0.08 g/m3) and (b) glaze ice (LWC ¼0.3 g/m3).
Fig. 12. Effect of icing time—AERTS S809 blade: (a) rime ice (LWC ¼0.05 g/m3) and (b) glaze ice (LWC ¼0.3 g/m3).
during AERTS case 14 is shown in Fig. 17(a) along with the baseline case of a clean blade. The corresponding increase in torque for this severe glaze ice condition was quantified to be
70.8%. AERTS cases 1 and 2 were tests under an extreme LWC condition (1.2 g/m3). Although this condition may be encountered very rarely by wind turbines in natural icing clouds, these tests serve as examples of worst-case scenarios. The recorded torque data can be of use for wind turbine blade design procedures as estimates for maximum loads and rotor imbalance. The corresponding torque rises for these two cases were measured to be 117.8% and 95.9% (see Fig. 17(b)). Tests were stopped when twice of the torque was reached due to potential shedding concerns.
66
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
4. Conclusions This work presents the first scaled ice accretion experiments on a rotating model wind turbine blade. The tip conditions of the test rotor were scaled to represent the 95% radial station of the NREL Phase VI Rotor under a variety of operating conditions. A modified Ruff method was used to scale the icing conditions. Tests were conducted in both the rime and glaze ice regimes. Rime ice shapes were obtained for icing times of up to 1.5 h (equivalent to 2.7 h full-scale) and compared very well with LEWICE predictions to within 1% of the accreted ice thickness. These results serve as a validation for the AERTS facility in reproducing representative icing clouds pertinent to wind turbines. Effects of icing parameter variations in AOA, temperature, LWC, and icing time were examined. The effect of AOA was small at lower values (21–41). The AOA effect on ice shapes obtained within the normal wind turbine blade operating range (41–81) was profound, which suggested that local AOA needs to be accurately predicted to ensure accurate prediction of ice accretion. While rime ice shapes were fairly independent of temperature, accreted glaze ice shapes showed a strong dependence on temperature. This was attributed to surface running water along the blade span for the glaze ice conditions due to the action of centrifugal forces. The thickness of rime ice shapes increased linearly with icing time, which stood in contrast to an irregular growth rate of glaze ice that was again attributed to surface running water. An increase in LWC resulted in a nearly linear increase in stagnation line ice thickness. These ice shapes corresponded to an icing time of up to 90 min in the rime ice regime. With respect to the glaze regime, the ice shape changes
Fig. 15. Torque measurement (Moderate Ice Coverage, glaze ice)—AERTS S809 blade.
disproportionately with icing time, a characteristic that is very difficult to predict with a 2D ice accretion model, as the 3D centrifugal force inherent in the rotating frame adds an additional difference to measured and predicted 2D ice shapes. Consequently, the measured glaze ice shapes did not compare well with LEWICE predictions. The total rotor torque was also monitored during the ice experiments, thus differentiating different classes of icing severity. The increases in torque requirement were in the orders of 5%, 25% and more than 70%, for the light-, moderate- and severe ice coverage, respectively. The tests conducted in the AERTS facility can serve as a database for future efforts in 3D ice accretion modeling by the wind energy community.
Fig. 17. Torque measurement (Severe Ice Coverage, glaze ice)—AERTS S809 blade: (a) LWC¼ 0.4 g/m3 and (b) LWC¼1.2 g/m3.
Fig. 16. Severe Ice Coverage (glaze ice).
Y. Han et al. / J. Wind Eng. Ind. Aerodyn. 109 (2012) 55–67
Acknowledgment This work was supported by startup funds for research in Wind Energy at the Pennsylvania State University. References Anderson, D., 2004. Manual of Scaling Methods, NASA/CR-2004-212875. Barber, S., et al., 2011. The impact of ice formation on wind turbine performance and aerodynamics. Journal of Solar Energy Engineering 133/011007, 1–9. Belte, D., 1987. In-flight ice accretion characteristics of rotor blade airfoil sections. In: AIAA-87-0176, AlAA 25th Aerospace Sciences Meeting. Reno, NV, 1987. Bourgault, Y., Boutanios, Z., Habashi, W.G., 2000. Three-dimensional eulerian approach to droplet impingement simulation using FENSAP-ICE, Part 1: model, algorithm, and validation. Journal of Aircraft 37 (1), 95–103. Byrkjedal, Ø., Vindteknikk, K., 2009. Estimating wind power production loss due to icing. In: Proceedings of the IWAIS XIII. Andermatt, 2009. Chattot, J.J., 2005. Helicoidal vortex model for steady and unsteady flows. Computers and Fluids 35 (7), 733–741. Finstad, K.J., Lozowski, E.P., Makkonen, L., 1988. On the median volume diameter approximation for droplet collision efficiency. Journal of the Atmospheric Science 45 (24), 4008–4012. Flemming, R.J., Lednicer, D.A., 1986. Correlation of icing relationships with airfoil and rotorcraft icing data. Journal of Aircraft 23 (10), 737–743. Fu, P., Farzaneh, M., 2010. A CFD approach for modeling the rime-ice accretion process on a horizontal-axis wind turbine. Journal of Wind Engineering and Industrial Aerodynamics 98, 181–188. ¨ ¨ Glauert, H., 1935. Das maximum der theoretisch moglichen ausnutzung des windes. Aerodynamic Theory (4), 169–360 (division L). Gray, V.H., von Glahn, U.H., 1958. Aerodynamic effects caused by icing of an unswept NACA 65A004 airfoil. National Advisory Committee for Aeronautics, NACA TN, Cleveland, OH 4155. Han, Y., Palacios, J., Smith, E., 2011. An experimental correlation between rotor test and wind tunnel ice shapes on NACA 0012 airfoils. In: Proceedings of the 2011-38-0092, SAE 2011 International Conference on Aircraft and Engine Icing and Ground Deicing. Chicago, 2011. Hand, M.M., Simms, D.A., Fingersh, L.J., Jager, D.W., Cotrell, J.R., Schreck, S., Larwood, S.M., 2001. Unsteady Aerodynamics Experiment Phase VI: Wind Tunnel Test Configurations and Available Data Campaigns. National Renewable Energy Laboratory, NREL/TP-500-29955. Golden, CO, 2001. Heinrich, A., et al., 1991. Aircraft Icing Handbook—Volume 1 of 3. Atlantic City International Airport, NJ 08405: Department of Transportation, Federal Aviation Administration Technical Center, DOT/FAA/CT-88/8-1. Hochart, C., Fortin, G., Perron, J., 2008. Wind turbine performance under icing conditions. Wind Energy (November), 319–333. Homola, M.C., Virk, M.S., Nicklasson, P.J., Sundsbø, P.A., 2012. Performance losses due to ice accretion for a 5 MW wind turbine. Wind Energy 15, 379–389.
67
Homola, M.C., et al., 2010. Effect of atmospheric temperature and droplet size variation on ice accretion of wind turbine blades. Journal of Wind Engineering and Industrial Aerodynamics 98 (12), 724–729 ) 2010. Homola, M.C., et al., 2011. Modeling of ice induced power losses and comparison with observations. In: Proceedings of the Winterwind 2011, International ˚ Sweden, 2011. Wind Energy Conference. Umea, Jasinski, W.J., Noe, S.C., Selig, M.S., Bragg, M.B., 1998. Wind turbine performance under icing conditions. Journal of Solar Energy Engineering 120, 60–65. Langmuir, I., Blodgett, K., 1946. A Mathematical Investigation of Water Droplet Trajectories. Army Air Forces Technical Report no. 5418. Makkonen, L., Laakso, T., Marjaniemi, M., Finstad, K.J., 2001. Modelling and prevention of ice accretion on wind turbines. Wind Engineering 25 (1), 3–21. Makkonen, L., Oleskiw, M.M., 1997. Small-scale experiments on rime icing. Cold Regions Science and Technology 25, 173–182. Manwell, J.F., McGowan, J.G., Rogers, A.L., 2010. Wind Energy Explained. John Wiley & Sons, New York, NY. Messinger, B., 1953. Equilibrium temperature of an unheated icing surface as a function of airspeed. Journal of the Aeronautical Sciences 20 (1), 29–42. Moriarty, P.J., Hansen, A.C., 2005. AeroDyn Theory Manual. Golden, CO. National Renewable Energy Laboratory. Nygaard, B.E.K., Kristja´nsson, J.E., Makkonen, L., 2011. Prediction of in-cloud icing conditions at ground level using the WRF model. Journal of Applied Meteorology and Climatology 50 (12), 2445–2459. Palacios, J., Han, Y., Brouwers, E., Smith, E., 2012. Icing environment rotor test stand liquid water content measurement procedures and ice shape correlation JAHS-1562. Journal of American Helicopter Society 57 (2) 022006-1-12. Reinert, T.F.R.J.N.R., Aubert, R.J., 2011. Oscillating airfoil icing tests in the NASA Glenn Research Center Icing Research Tunnel. In: Proceedings of the SAE 2011 International Conference on Aircraft and Engine Icing and Ground Deicing. Chicago, IL, 2011. Ruff, G.A., 1986. Analysis and Verification of the Icing Scaling Equations. AEDC-TR85-30, vol. 1. (Rev.). Ruff, G.A., Berkowitz, B.M., 1990. Users Manual for the NASA Lewis Ice Accretion Prediction Code (LEWICE), NASA CR 185129. Schmitz, S., 2011. XTurb-PSU—A Wind Turbine Design and Analysis Tool. The Pennsylvania State University, University Park. Shin, J., Berkowitz, B., Chen, H.H., Cebeci, T., 1994. Prediction of ice shapes and their effect on airfoil drag. Journal of Aircraft 31 (2), 263–270. Somers, D.M., 1997. Design and Experimental Results for the S809 Airfoil. National Renewable Energy Laboratory, Golden, CO. Szasz, R.Z., Fuchs, L., 2011. Ice accretion prediction on wind turbines based on a combined LES-LPT method. In: Proceedings of the Winterwind 2011, Interna˚ Sweden, 2011. tional Wind Energy Conference. Umea, Tran, P., et al., 1994. Ice accretion on aircraft wings with thermodynamic effects. Journal of Aircraft 32 (2), 444–446. Wright, W.B., 2005. Validation results for LEWICE 3.0. In: Proceedings of the AIAA2005-1243, 43rd Aerospace Sciences Meeting and Exhibit. Reno, Nevada, 2005. Wright, W.B., 2008. User’s Manual for LEWICE Version 3.2, NASA CR 2008-214255.