Dynamic Stall Model for Tower Shadow Effects on Downwind Turbines and Its Scale Effects

Abstract

A dynamic stall model for tower shadow effects is developed for downwind turbines. Although Munduate’s model shows good agreement with a 1.0 m wind tunnel test model, two problems exist: (1) it does not express load increase before the entrance of the tower wake, and (2) it uses the empirical tower wake model to determine the wind speed profile behind the tower. The present research solves these problems by combining Moriarty’s tower wake model and the entrance condition of the tower wake. Moriarty’s model does not require any empirical parameter other than tower drag coefficient and it expresses positive wind speed around the tower also. Positive wind speed change is also allowed as the tower wake entrance condition in addition to the negative change observed in the previous model. It demonstrates better agreement with a wind tunnel test and contributes to the accuracy of the fatigue load, as it expresses a slight increase in load around the entrance of the tower wake. Furthermore, the scale effects are also evaluated; lift deviation becomes smaller as the scale increases, i.e., lower rotor speed.

1. Introduction

Upwind turbines, which have rotors in front of the towers, have been predominant throughout the decades of commercial wind turbine history. However, downwind turbines, which have rotors behind the towers, are gathering attention due to their technical and economic advantages regarding modern/future wind turbines. Downwind rotors generally have negative rotor tilt to avoid collision of the blade and the tower. One of the advantages of downwind turbines is performance in complex terrains. The negatively tilted rotors of downwind turbines are favorable for upflow wind in complex terrain [1]. Furthermore, the yaw measurement devices in front of the rotor and the nacelle also contribute to yaw control accuracy [2]. Downwind turbines can be advantageous in floating offshore wind turbine (FOWT) applications, so the share of the downwind rotor is much higher in FOWTs [3,4,5] compared to the rest of the market. New FOWT concepts are now appearing, with some downwind rotors [6,7] inclining rearward a couple of degrees due to the rotor thrust, causing the rotor–wind misalignment to become smaller in complex terrain conditions. Rotor position and rotor coning further contribute to the stability of floating turbines with small yaw stiffness and damping [8,9]. Downwind turbines are also considered to be advantageous for large-scale wind turbines, mainly due to their compatibility for lighter, more elastic blades and load mitigation by appropriately coned rotors [10].

The most essential drawback of downwind turbines is the tower shadow effect, which generates impulsive loads and infrasound when the blades pass through the wake of the tower [11]. Many design load cases are defined by combinations of wind turbine conditions, as well as various wind and marine conditions based on international design standards, such as IEC61400-1 [12] and IEC61400-3-1 [13]. All of the aerodynamics, the elasticities of the structure, the controls, and the hydrodynamics strongly affect the loads and performances of the wind turbines. The blade-element momentum (BEM) method is commonly used in analyses, with some modifications or extensions for three-dimensional and dynamic effects [14,15], as it is accurate enough with low cost analysis. Therefore, development of the tower shadow model for BEM is the most important technical challenge in the design and analysis of downwind wind turbines. Previous studies considered the variable loads of downwind turbines by tower shadow effects, demonstrating applications of the tower wake wind speed profile while ignoring the interaction between the rotor and the tower. Matiz-Chicacausa and Lopez [16] conducted the analysis of the tower shadow effects by the actuator line model, which showed good agreement with computational fluid dynamics (CFD). Wang and Coton [17] developed a high-resolution tower shadow model, which showed good agreement with an experiment except for high angle of attack conditions. Zahle et al. [18] conducted a two-dimensional CFD experiment for tower shadow effects on three different tower configurations. Zhao et al. [19] compared upwind and downwind rotors with two different rotor speed conditions using CFD. Van der Male et al. [20] showed that the tower shadow of downwind turbine strongly affects the fatigue damages of downwind turbine by the aeroelastic simulation with Madsen’s wake profile model. The aerodynamic interaction with between the rotor and the tower and dynamic effects are ignored in the study. Yoshida and Kiyoki [21] developed the load equivalent tower shadow modeling method for the BEM of downwind turbines. It defines a bell-shaped wind speed profile of the tower wake by its three parameters, i.e., depth and width of the tower wake profile and the defined point, which could be adapted to the load history through the tower wake using the wind turbine CFD. This was a considerable finding, as it considered the aerodynamic interaction between the rotor and the tower, which could not be taken into account in previous models, thereby providing realistic load fluctuation in the tower shadow. This technique was practically used in the design and analysis of commercial downwind turbines, such as SUBARU 80/20 [22], later Hitachi 2 MW, and Hitachi 5 MW [23]. This method is still useful for analysis but not for practical design applications, as it needs CFD for each condition to identify the parameters. In addition, this model does not consider the dynamic tower shadow effects. Munduate et al. [24] developed a dynamic tower shadow model, which considers dynamic stall effect. Although it ignores the mutual interaction between the rotor and the tower, it shows good agreement with wind tunnel tests with a 1.0 m rotor model, particularly in the context of asymmetry between the entrance and exit of the tower wake. However, the model still demonstrated two problems: (1) it did not express the load increase before the entrance of the tower wake, thereby affecting the fatigue, and (2) it uses an empirical tower wake model to determine the wind speed profile behind the tower.

Considering these situations, a dynamic tower shadow modeling method is developed herein for the BEM calculation of downwind turbines. Munduate’s model is modified and extended to solve problems (1) and (2), as mentioned above. Furthermore, the scale effect of the model is also discussed in this study.

2. Methodology

This research is an extension of Munduate et al. [24], aiming for better agreement at the entrance of the wake without using any empirical parameters. The method consists of the following steps:

(1)

The dynamic tower shadow model;

(2)

The wind speed profile model behind the tower;

(3)

The threshold of the entrance of the tower wake.

There are some assumptions in the present methodology. The tower wake is represented by the average wind speed profile with the turbulence ignored, and the trailing vorticity does not affect the lift of the blade element.

2.1. Dynamic Tower Shadow Model

The dynamic tower shadow model is based on Munduate et al. [24]. The deviation of the lift coefficient of the blade section ΔC_l while the blade section is passing through the tower wake is

$$\mathrm{\Delta }{C}_l=\frac{2\pi }{W}\left[w_g\left(0\right)\mathsf{\Psi }\left(s\right)+{{\displaystyle \int }}_0^s\frac{d{w}_g}{ds}\mathsf{\Psi }\left(s-\sigma \right)d\sigma \right]$$

(1)

where W is the relative wind speed of the blade section without the tower effect, and w_g is the wind speed profile of the rotor plane without the rotor, which is discussed in the next section. In addition, s is the normalized time defined as below according to the blade chord length c and time from the entrance of the tower wake t.

$$s=\frac{2Wt}{c}$$

(2)

The definition of the entrance of the tower wake is discussed in Section 2.3.

Sears and Sparkes [25] proposed the following equation for the Kussner function Ψ:

$$\mathsf{\Psi }\left(s\right)=1-A_{1}exp\left(-b_{1}s\right)-A_{2}exp\left(-b_{2}s\right)$$

(3)

$$A_1=A_2=0.5, b_1=0.13, b_2=1.0$$

Therefore, Equation (1) is

$$\mathrm{\Delta }{C}_l=\frac{2\pi }{W}\left[w_g\left(s\right)-X\left(s\right)-Y\left(s\right)\right]$$

(4)

$$\left\{\begin{array}{c}X\left(s\right)=A_1{{\displaystyle \int }}_0^s\frac{d{w}_g}{ds}exp\left[-b_1\left(s-\sigma \right)\right]d\sigma \\ Y\left(s\right)=A_2{{\displaystyle \int }}_0^s\frac{d{w}_g}{ds}exp\left[-b_2\left(s-\sigma \right)\right]d\sigma \end{array}\right.$$

2.2. Wind Speed Profile behind the Tower

Tower wake models provide wind speed profiles behind the tower. Many models have previously been proposed, as shown in

Table 1. Tower wake models.

Model	Outlines	Equation
Blevin [26]	$\frac{U}{U_0}=1-cexp\left(-0.69\frac{y^2}{b^2}\right)$	(5)
Madsen [27]	$\frac{U}{U_0}=\frac{\sqrt{3}}{2}\sqrt{\frac{K\sigma }{x}}\left[1-tan{h}^2\left(\eta \right)\right]$	(6)
Moriarty [28]	$\frac{U}{U_0}=\frac{U_{Ex}}{U_0}-\frac{C_d}{\sqrt{\delta }}{cos}^2\left(\frac{\pi }{2}\frac{\eta }{\sqrt{\delta }}\right)$	(7)
Munduate [24]	$\frac{U}{U_0}=-\frac{D_v}{2}\left[1-cos\frac{2\pi \left({\Psi }_t-\frac{\pi }{2}\right)}{2{\Psi }_0}\right]$	(8)
Powles [29]	$\frac{U}{U_0}=1-\Delta {cos}^2\left(\frac{y}{wD}\pi \right)$	(9)
Schlichting [30]	$\frac{U}{U_0}=\frac{C_d}{4\sqrt{\pi }}\sqrt{\frac{U_{0}D}{\nu }}{\left(\frac{x}{D}\right)}^{-\frac{1}{2}}exp\left(-\frac{y^2{U}_0}{4x\nu }\right)$	(10)

.

Although explanations of each parameter are omitted here, it is common for all of the models above that they do not consider the influence of the rotor. These models are categorized into several groups. Powles’ model includes a form of a square of the cosine. Munduate’s and Blevin’s models show similar forms, which contain one or more empirical parameters. Schlichting’s also has a similar form but is more useful as it does not contain any empirical parameter. Madsen’s model shows another form as the square of hyperbolic tangent, including an empirical parameter. Moriarty’s model, which is used in the aeroelastic simulation code FAST/Aerodyn [15], is applied in this research. However, it is also an extension of Powles’ model; it expresses the local wind speed increase around the tower and does not contain any empirical parameter, unlike other models.

The formulation of Moriarty’s model is as below. The wind speed profiles outside (U_Ex) and inside of the wake (U) at the longitudinal and lateral positions (x, y) behind the tower center are provided. Here, the wind speeds are normalized by the free stream wind speed U₀ and the tower radius R, respectively. C_d represents the drag coefficient of the tower section.

$$\left\{\begin{array}{cc}\frac{U_{Ex}}{U_0}& =1-\frac{{\xi }_C{}^2-{\eta }^2}{{\left({\xi }_C{}^2+{\eta }^2\right)}^2}+\frac{C_d}{2\pi }\frac{{\xi }_C}{{\xi }_C{}^2+{\eta }^2}\cdots \left|\eta \right|>\sqrt{\delta }\\ & \frac{U}{U_0}=\frac{U_{Ex}}{U_0}-\frac{C_d}{\sqrt{\delta }}{cos}^2\left(\frac{\pi }{2}\frac{\eta }{\sqrt{\delta }}\right)\cdots \left|\eta \right|\le \sqrt{\delta }\end{array}\right.$$

(11)

where

$$\begin{array}{c}\delta =\sqrt{{\xi }^2+{\eta }^2}\\ \xi =x/R\\ \eta =y/R\\ {\xi }_C=\xi +0.01\end{array}$$

2.3. Threshold of the Entrance of the Tower Shadow

The entrance of the tower shadow is defined in the present dynamic tower shadow model in Section 2.1. More practically, it is defined as the initial point of s in Equations (1) and (4), according to Munduate et al. [24], i.e.,

$$\frac{U}{U_0}<1$$

(12)

The model does not show the load increase around the entrance of the tower wake as mentioned in Section 1, as its wind speed profile does not express the wind speed increase around the tower, as seen in Equation (8). However, the present model can show the wind speed increase around the entrance of the tower wake as Moriarty’s model is applied. The following condition is proposed for the tower wake entrance condition in this research, to deal with both the positive and negative deviation with the 1% of dead band:

$$\left|\frac{U}{U_0}-1\right|>0.01$$

(13)

The sensitivity to the width of the dead band is discussed in Section 3.6.

3. Validation of the Model

The model in the previous section is validated herein by the wind tunnel test conducted by Munduate et al. [24].

3.1. Model

The wind tunnel test model and its installation in the wind tunnel are shown in Figure 1 and

Table 2. Specifications of the wind tunnel test model described in [24].

Item	Specification
Rotor diameter	1.0 m
Number of blades	2
Blade section	NACA4415
Chord length c	0.1 m
Twist angle	0 deg
Blade set angle	12 deg
Rotor tilt angle	0 deg
Rotor coning angle	0 deg
Dummy tower diameter D	0.07 m
Dummy tower position	0.14 m

. The test was conducted in the Low Speed Wind Tunnel at the University of Glasgow (UG), which has a test section with a 2.13 m width and a 1.61 m height. A total of 23% of the blockage factor can affect the free stream wind speed up to 6–12% in cases of 0.6–0.9 of the thrust coefficients [31]. However, no correction is applied for the influence in this study, as it affects only the operation condition in the discussion of the variation of the lift coefficients. The variation of the lift coefficient measured by the pressure transducers at 75% rotor radius was used in the experiment.

3.2. Test Conditions

The wind tunnel test was conducted in the three cases as shown in

Table 3. Test conditions [24].

Case	1	2	3
Wind speed U0	9.0 m/s	11.0 m/s	11.7 m/s
Rotor speed n	918 min−1	1122 min−1	1188 min−1
Tip speed ratio λ	5.34	5.33	5.31
Tangential speed at ηB = 75% VT	36.0 m/s	44.1 m/s	46.7 m/s

. The tip speed ratios were maintained at about 5.3 in three different wind speeds.

3.3. Wind Speed Distribution

The wind speed distribution calculated in this study is shown in Figure 2. Here, the drag coefficient of the tower section C_d is assumed to be 1.2, in accordance with the Reynolds number (4.5 × 10⁴) [32]. It shows narrow and deep distribution, and the local wind speed is higher than the free stream just outside of the tower wake at x/R = 2.0. It becomes shallower and wider as the wind flows downstream. Here, the rotor position of the test model is at x/R = 4.0.

3.4. Validation of the Model

The wind speed distribution on the rotor plane normalized by free stream wind speed to rotor azimuth and the blade station radius is shown in Figure 3. Here the influence of the rotor is not considered. The normalized wind speeds are identical in both figures, showing symmetrical distribution around 180°, where the first blade is located just behind the tower. Moreover, it is geometrically clear that the inner sections are affected by the tower shadow in wider azimuth position.

Distributions of the deviation of the lift coefficients calculated by the present method by Munduate’s and the present thresholds are shown in Figure 4 and Figure 5, for the three cases in Table 3. Although the overall characteristics such as the depth and the lift drop are almost same, they are slightly different due to differences in the entrance conditions around 160°–180°.

The deviations of the lift coefficient calculated by the present method with Munduate’s threshold are shown in Figure 6. This shows the wind tunnel test (“Experiment UG”), static (“Steady UG”) and dynamic (“Unsteady UG”) simulation results as seen in [24], as well as the present static (“Steady KU”) and dynamic (“Unsteady KU”) simulation results. Here, UG and KU stand for University of Glasgow and Kyushu University. The steady data are calculated by X = Y = 0 using Equation (4). “Unsteady KU” is identical to the data at η_B = 75%, as shown in Figure 4. Both unsteady methods show much better agreement with the test data compared to the steady data. Although the deviation caused by the differences in the wind speed profiles is small, better agreement is demonstrated by the present method with Munduate’s threshold. However, neither of these trials show an increase in lift before it drops at around 160°–180°. In the same manner, the deviations of the lift coefficient calculated in accordance the present method with the present threshold are shown in Figure 7. Unlike Munduate’s threshold, it shows a slight increase in lift at the entrance of the tower wake, as does the experiment. However, it is less important than the range from the fatigue and the extreme loads points of view, the recoveries from the lift drop are slower than the experiment. This is considered as the future problem in this research.

The maximum and the minimum of the three cases are shown in Figure 8. The maximum values by Munduate are zero, as shown above, whereas the present method shows slightly positive values, as does the experiment. The minimum values also show better agreement with the experiment. The ranges between the maximum and the minimum are shown in Figure 9. The present method shows better agreement with the experiment as compared with Munduate’s. This difference is quite important, as the range strongly affects the accuracy of the fatigue load.

3.5. Sensitivity of Drag Coefficient of the Tower Section

One of the advantages of the present method is that no empirical parameter is required. The tower drag coefficient is the only parameter explicitly assumed. The tower drag coefficient varies in accordance with the Reynolds number, as well as the change of the cross sections and the surface conditions. To evaluate the sensitivity of the present model to the tower drag coefficient, the simulation results using C_d = 0.6 (−50 %) and 1.8 (+50 %) at V_T = 36.0 m/s are shown in Figure 10 and Figure 11. The tower drag coefficient strongly affects the depth of the lift drop, as well as the increase in the lift around the entrance of the wake. Therefore, the accuracy of the tower drag coefficient is quite important in the present model.

3.6. Sensitivity of the Dead Band of Wake Entrance Condition

The wake entrance condition in the present research is shown in Equation (13). The sensitivity of the dead band is studied here. Those calculated by 0.5% (Figure 12) and 1% (Figure 7) of the dead band are almost identical. However, in cases of 2% of the dead band (Figure 13), the increase in lift are a little steeper than those of 1%. Considering the results, the dead band is determined as 1 % from Equation (13) in this study.

4. Scale Effect

The present method shown in Section 2 is based on the similarity of the normalized time s shown in Equation (2). The chord length c is proportional to the scale of the rotor, whereas the inflow wind speed W is almost same at the same wind speed and the tip speed ratio. The scale effect of s is evaluated in this section.

4.1. Analysis Outline

The wind tunnel test model shown in the previous chapter is herein scaled up to 300% and 1000%. Specifications of the three configurations are summarized in

Table 4. Specifications of the analysis of scale effects.

Scale	100%	300%	1000%
Rotor diameter	1.0 m	3.0 m	10.0 m
Number of blades	2	2	2
Blade section	NACA4415	NACA4415	NACA4415
Chord length	0.1 m	0.3 m	1.0 m
Twist angle	0 deg	0 deg	0 deg
Blade set angle	12 deg	12 deg	12 deg
Rotor tilt angle	0 deg	0 deg	0 deg
Rotor coning angle	0 deg	0 deg	0 deg
Tower diameter	0.07 m	0.21 m	0.7 m
Tower position	0.14 m	0.42 m	1.4 m
Reynolds number of the tower	4.5×104	1.4 ×105	4.5×105
Drag coefficient of the tower section [32]	1.2	1.2	0.33 (1.2)

. The 100% model is identical to the model presented in Section 3.2. Appropriate dimensions of blade, tower, and their relative position are not discussed in this study, assuming their similar configuration. In addition, the Reynold number effect on the blade is ignored, as this study focuses on the effect of the blade passing speed.

4.2. Analysis Conditions

Specifications of the models are shown in Table 4. The drag coefficients of the tower are estimated by reference [32]. C_d = 1.2 is also analyzed for the 1000% model for the comparison with other models. The analysis conditions are shown below.

Wind speed U₀: 9.0 m/s.

Tip speed ratio λ: 5.3.

4.3. Analysis Results

Distributions of wind speed throughout the rotor plane and the deviation of the lift coefficient are shown in Figure 14 and Figure 15a for the 300% and the 1000% models at C_d = 1.2, respectively. Figure 5a shows the same conditions for the 100% model. The deviations of the lift coefficients are completely different, demonstrating that the tower shadow effect becomes much smaller as the blade passing speed increases. The increase in the load around the entrance of the tower shadow also becomes minor.

The 1000% model with the Reynolds number effect on the tower drag is shown in Figure 12b. The drag of the cylinder in this Reynolds number, just above the critical Reynolds number, is quite small—C_d = 0.33. The deviations of the lift are also fairly smaller as C_d is remarkably smaller.

Although these results are dependent on the configurations and operational conditions, the dynamic tower shadow effect generally becomes smaller as the rotor speed decreases.

5. Conclusions

Although Munduate’s dynamic tower shadow model shows good agreement with the experiment using the 1.0 m rotor, two problems persist: (1) it does not express the load increase before the entrance of the tower wake, thereby affecting the fatigue, and (2) it uses the empirical tower wake model to determine the wind speed profile behind the tower.

The present research solves these problems by applying Moriarty’s tower wake model and the tower wake entrance conditions obtained from negative wind speed to accept both negative and positive deviations, demonstrating better agreement with the experiment. The present model may be effective in improving the accuracy of fatigue load analysis, as the range between the maximum and minimum load fluctuations around a tower is vital. The results show that load fluctuation decreases as the scale is increased.