Numerical Analysis of the Effects of Rotating Wind Turbine Blades on the Aerodynamic Forces Acting on Tower

We have investigated the effects of the rotating blades of an upwind-type three-blade horizontal-axis wind turbine (HAWT) on the basic characteristics of aerodynamic forces acting on its tower by conducting improved delayed detached-eddy simulations (DESs). Three tip-speed ratios were considered for the operating conditions of the HAWT: λ = 3 (low), λ = 6 (optimum), and λ = 10 (high). The diversion of the flow approaching the tower by the rotating blades and the low-pressure region that formed downwind of the blades significantly affected the aerodynamic forces acting on the tower. For example, the azimuth angle around the tower at which the pressure reached a maximum at each height shifted significantly in the direction of the movement of the blade passing the tower because of the diversion of the flow by the blades. Fluctuations in the lift force of the tower were significantly larger than those in its drag force because of the low-pressure region downwind of the blades.


Introduction
When designing a wind turbine, it is important to avoid the excitation of resonant oscillations of the wind turbine tower by rotor thrust fluctuations at the blade-passing frequency or, to a lesser extent, at the rotational frequency.Wind turbine towers are customarily divided into three categories according to their stiffness: stiff, soft, and soft-soft [1][2][3].A stiff tower is one whose fundamental natural frequency is above the blade-passing frequency range, a soft tower is one whose natural frequency is between the blade-passing frequency range and the rotor frequency range, and a soft-soft tower is one whose natural frequency is below both the rotor frequency range and the blade-passing frequency range.Because stiff towers usually require extra material not otherwise required for strength, soft towers are generally preferred.Soft and soft-soft towers can however be excited during the startup or shutdown of the turbine.In addition, dynamic magnification directly impacts the fatigue loads.Therefore, it is essential to design wind turbines and towers with the proper fatigue strength by considering their vibratory characteristics.
One of the main causes of rotor-thrust fluctuations at the blade-passing frequency is blade-tower interaction (BTI).Regardless of whether the configuration is upwind or downwind, the thrust force on the rotor of a horizontal-axis wind turbine (HAWT) decreases when a blade passes by its tower Energies 2017, 10, 121 2 of 18 because the wind velocity of the wind approaching the tower and the wind velocity in the wake of the tower are reduced owing to the existence of the tower.The effects of the BTI on the rotor thrust have been extensively studied with computational fluid dynamics (CFD) simulations for upwind-type HAWTs [4][5][6][7][8][9][10][11] and downwind-type HAWTs [12,13].However, studies on the effects of the BTI on the aerodynamic forces acting on the tower are very limited [11].These studies involved CFD simulations of the National Renewable Energy Laboratory (NREL) Phase IV turbine [14] using Reynolds-averaged Navier-Stokes (RANS) turbulence models.Their main findings are as follows: Gómez-Iradi et al. [5] showed that the pressure on the upwind surface of the tower was reduced above the lowest point of the rotor at a certain rotational speed.This effect was evident even when the tips of the two blades were the maximum distance away from the tower (i.e., the angle between the tower and the blade was 90 • ).Wang et al. [9] reported that the pressure drop on the upwind surface of the tower, which Gómez-Iradi et al. [5] showed, was very weak at a certain tip-speed ratio when the tips of the two blades were relatively far from the tower.Using the HAWT with a downwind configuration, Zahle et al. [13] showed that the rotor had a strong effect on the tower shedding frequency, causing vortex lock-in at the upper part of the tower under certain flow conditions.They suggested the possibility that these lock-in phenomena were due to numerical artifacts of the RANS turbulence model.They also emphasized the necessity of further investigation using a large-eddy simulation (LES) or a hybrid RANS-LES model, which can reproduce separated flows more realistically.
In this study, as a first step toward clarifying the effects of aerodynamic forces acting on the tower of a wind turbine on the vibratory characteristics of the tower, we investigated the effects of the rotating blades of an upwind-type HAWT on the basic aerodynamic characteristics of the tower (e.g., the drag, lift, and pressure coefficients).This was done by conducting CFD simulations with the improved delayed detached-eddy simulation (IDDES) model [14], which is a hybrid RANS-LES model, of the wind flow around the wind turbine and its tower at low, optimum, and high tip-speed ratios.

Numerical Approach
In this study, a right-handed Cartesian coordinate system (x 1 , x 2 , x 3 ) = (x, y, z) is employed, in which the z-direction is aligned in the vertical direction.The CFD software utilized to simulate the wind flow field is ANSYS Fluent 17.1 (ANSYS, Inc., Canonsburg, PA, USA) [15,16].

Wind Turbine and Tower
Figure 1 shows the geometry of the wind turbine and tower investigated by the CFD simulation.The rotor and nacelle are modeled as the upwind-type HAWT used by Eriksen and Krogstad [17] and Krogstad and Eriksen [18] in their wind tunnel experiments.The diameter of the rotor is D = 0.894 m, and the center of the rotor is located 0.817 m above the floor.The blades have a non-uniform twist and chord-length distribution, and have an NREL S826 airfoil profile [19].The tower has a uniform diameter of d = 0.061 m in the vertical direction.It is divided into 10 segments, and the aerodynamic forces acting on each segment are analyzed.The segments have identical lengths of h = 0.0772 m.

Governing Equations and Discretization Method
The detached-eddy simulation (DES) model [20] is a hybrid RANS-LES model that treats near-wall regions in a RANS-like manner, and treats the rest of the flow in an LES-like manner.This model has the potential to save a large amount of computing resources compared with a pure LES model.However, modeled-stress depletion and grid-induced separation are problems in this model.The IDDES model mitigates these effects by delaying the RANS-to-LES transition and increasing the modeled stress contribution across the interface.The flow field around the wind turbine is assumed to be incompressible and isothermal in air at a temperature of 15 °C.The governing equations for the IDDES model are the continuity equation: the Navier-Stokes equation: the transport equation for the turbulence kinetic energy k: and the transport equation for the specific dissipation rate ω: ( ) where ui is the wind-velocity component in the xi direction; p is the pressure; ν is the kinematic viscosity; μt is the turbulence viscosity; t is the time; ρ is the air density; Sij is the mean strain rate; LIDDES is the IDDES length scale [14]; F1 is a blending function [21]; σk and σω are the turbulent Prandtl numbers for k and ω, respectively; and β is a model constant [22].μt is computed as follows: The flow field around the wind turbine is assumed to be incompressible and isothermal in air at a temperature of 15 • C. The governing equations for the IDDES model are the continuity equation: the Navier-Stokes equation: the transport equation for the turbulence kinetic energy k: and the transport equation for the specific dissipation rate ω: where u i is the wind-velocity component in the x i direction; p is the pressure; ν is the kinematic viscosity; µ t is the turbulence viscosity; t is the time; ρ is the air density; S ij is the mean strain rate; L IDDES is the IDDES length scale [14]; F 1 is a blending function [21]; σ k and σ ω are the turbulent Prandtl numbers for k and ω, respectively; and β is a model constant [22].µ t is computed as follows: Energies 2017, 10, 121 where α* is a low-Reynolds number correction coefficient [21], F 2 is a blending function [21], and α 1 is a model constant.
The governing equations are discretized by the finite-volume method.The advection terms of the Navier-Stokes equations are discretized by the bounded central-difference scheme.The advection terms of the transportation equations for k and ω are discretized by the QUICK scheme.Other spatial derivatives are discretized by the central-difference scheme.The time integration is performed using the first-order implicit method.

Numerical Setup
The computational grids, the boundary conditions, and the total number of computational cells for all run cases are shown in Figure 2 and Table 1.The origin of the coordinate system is defined as the intersection point between the floor and the vertical line intersecting the center of the rotor.The computational domain consists of the rotational domain, which includes the rotor, and the non-rotating exterior domain.The size of the non-rotating exterior domain is 12.5D × 3D × 2D, which is approximately the same as that of the wind tunnel used by Krogstad and Eriksen in their experiments [17,18].The diameters of the rotational domain D R for all run cases are listed in Table 1.
Here, d tr (=(1.13D− D)/2) is the shortest distance between the tip of a blade and the edge of the rotational domain.The blades, tower, and wind-tunnel walls are covered with boundary-layer meshes.Ideally, the first grid nodes over the surface of the blades should be set as y+ < 1 to resolve the eddies in the semi-viscous near-wall region.However, this requires an extremely large number of computational cells, in particular at a high tip-speed ratio.where α* is a low-Reynolds number correction coefficient [21], F2 is a blending function [21], and α1 is a model constant.
The governing equations are discretized by the finite-volume method.The advection terms of the Navier-Stokes equations are discretized by the bounded central-difference scheme.The advection terms of the transportation equations for k and ω are discretized by the QUICK scheme.Other spatial derivatives are discretized by the central-difference scheme.The time integration is performed using the first-order implicit method.

Numerical Setup
The computational grids, the boundary conditions, and the total number of computational cells for all run cases are shown in Figure 2 and Table 1.The origin of the coordinate system is defined as the intersection point between the floor and the vertical line intersecting the center of the rotor.The computational domain consists of the rotational domain, which includes the rotor, and the non-rotating exterior domain.The size of the non-rotating exterior domain is 12.5D × 3D × 2D, which is approximately the same as that of the wind tunnel used by Krogstad and Eriksen in their experiments [17,18].The diameters of the rotational domain DR for all run cases are listed in Table 1.
Here, dtr (=(1.13D− D)/2) is the shortest distance between the tip of a blade and the edge of the rotational domain.The blades, tower, and wind-tunnel walls are covered with boundary-layer meshes.Ideally, the first grid nodes over the surface of the blades should be set as y+ < 1 to resolve the eddies in the semi-viscous near-wall region.However, this requires an extremely large number of computational cells, in particular at a high tip-speed ratio.According to grid-sensitivity studies on the NREL S826 airfoil, Krogstad and Lund [23] set the first grid nodes over the surface of the blades to be y+ < 5 in their CFD simulations for the same HAWT rotor that is used in the present study.They utilized the sear stress transfer (SST) k-ω turbulence model with low-Reynolds number corrections [21], which is the same turbulence model used in this study in the near-wall region.Because their CFD results for the power coefficients and thrust coefficients well matched their experimental results over a wide range of tip-speed ratios, we also set the first grid nodes over the surface of the blades to be y+ < 5 in all run cases by setting them as 2 × 10 −5 m off the surface of the blades.
The first grid nodes over the surface of the tower are also set at 2 × 10 −5 m off the surface of the tower, resulting in y+ < 1.As shown in Figure 2c, structured meshes are used in a region around and downwind of the tower.This region extends from z = 0 to z = 9.9 h.The numbers of the cells in this region are listed in Table 1 for all run cases.Here, case G_fine has a number of cells that is approximately 1.5 times larger in each direction (x, y, and z) for the rectangular mesh and in each radius, θ, and z for the O-type mesh than case λ_6.In addition, case λ_6 has a number of cells that is approximately 1.5 times larger each direction (x, y, and z) for the rectangular mesh, and in each radius, θ, and z for the O-type mesh than case G_coarse.
At the inlet boundary, a stream-wise wind velocity of U ref = 10 m/s with a turbulence intensity of TI = 0.3% is implemented.At the outlet boundary, the outflow boundary condition is imposed.On the surfaces of the blades and tower and on the wind-tunnel walls, the no-slip boundary conditions are set.The sliding mesh technique is used to couple the rotational domain and the stationary domain.The direction of the rotor rotation is counterclockwise when viewed from upwind (Figure 1b).By changing the rotational speed of the rotor Ω, the tip-speed ratio λ (=DΩ/2U ref ) is set at 3, 6, or 10.Except for cases dt_fine and dt_coarse, the time-step sizes are calculated as dt = 2 • /Ω.This yields time-step values of 5.20 × 10 −4 s when λ = 3, 2.60 × 10 −4 s when λ = 6, and 1.56 × 10 −4 s when λ = 10.The time-step sizes for cases dt_fine and dt_coarse in Table 1 are calculated as dt = 1 • /Ω and dt = 4 • /Ω, respectively.The statistics of the flow are calculated over the time interval t ≈ 0.9-2.8s, which is approximately 328 in non-dimensional time based on d.This provides a sufficient interval for evaluating the flow statistics around a circular cylinder [24].The time-averaged value of an arbitrary physical quantity ϕ is indicated as ϕ. ) evaluated from the simulations to those obtained from the wind-tunnel experiments [18].Here, Q is the torque on the rotor, A is the rotor swept area, and T is the total thrust

Validation of Numerical Approach
Energies 2017, 10, 121 6 of 18 on the blades.The simulations provide a reasonable qualitative and quantitative approximation of the experimentally observed results for both C P and C T .The discrepancies of C P and C T between case λ_6 and other cases with λ = 6 are very small and are at most 0.0067 (dt_coarse) and 0.0048 (G_coarse), respectively.
Energies 2017, 10, 121 6 of 18 and CT between case λ_6 and other cases with λ = 6 are very small and are at most 0.0067 (dt_coarse) and 0.0048 (G_coarse), respectively.) of each tower segment.Here, Fx and Fy are the drag and lift forces, respectively, acting on each tower segment, and σFx and σFy denote the root mean squares of Fx and Fy, respectively.Except for case G_coarse, which is not shown in the figures, the discrepancies of these drag and lift coefficients between case λ_6 and other cases with λ = 6 are small.and CT between case λ_6 and other cases with λ = 6 are very small and are at most 0.0067 (dt_coarse) and 0.0048 (G_coarse), respectively.) of each tower segment.Here, Fx and Fy are the drag and lift forces, respectively, acting on each tower segment, and σFx and σFy denote the root mean squares of Fx and Fy, respectively.Except for case G_coarse, which is not shown in the figures, the discrepancies of these drag and lift coefficients between case λ_6 and other cases with λ = 6 are small.Figure 6 shows the sensitivity of D R , the grid resolution of the structured mesh region near the tower, and dt on the horizontal distributions of the mean pressure coefficient and the fluctuating pressure coefficient C p ' (=2σ p /ρU ref 2 ) on the tower surface at the mid-height of Segment 6.Here, p ref is the mean static pressure at (−4.1D, 0, 0.91D), σ p denotes the root mean square of p, and θ is the azimuth angle around the tower, whose direction is defined in Figure 1d.Except for case G_coarse, which is not shown in the figures, the discrepancies of these pressure coefficients between case λ_6 and other cases with λ = 6 are small.This tendency is observed at other heights.The aforementioned comparisons with the wind-tunnel experiments and the aforementioned sensitivity analyses indicate that the grid resolution, the rotational-domain size, and the time-step sizes for cases λ_3, λ_6, and λ_10 are reasonable for discussing the characteristics of the aerodynamic forces acting on the tower.In the following sections, the results for cases λ_3, λ_6, and λ_10 are discussed.

Pressure Coefficients on the Tower
Figure 7 shows the distributions of C p on the tower surface at the mid-heights of each tower segment.When λ = 10, the distributions of C p around the tower are almost symmetric with respect to θ = 0 for all segments.However, when λ = 3 and λ = 6, the angle θ at which C p is maximized shifts in the −θ direction at upper segments, particularly at segments 1 and 2. This shifting is caused by the diversion of the flow approaching the tower by the rotating blades.) on the tower surface at the mid-height of Segment 6.Here, pref is the mean static pressure at (−4.1D, 0, 0.91D), σp denotes the root mean square of p, and θ is the azimuth angle around the tower, whose direction is defined in Figure 1d.Except for case G_coarse, which is not shown in the figures, the discrepancies of these pressure coefficients between case λ_6 and other cases with λ = 6 are small.This tendency is observed at other heights.
The aforementioned comparisons with the wind-tunnel experiments and the aforementioned sensitivity analyses indicate that the grid resolution, the rotational-domain size, and the time-step sizes for cases λ_3, λ_6, and λ_10 are reasonable for discussing the characteristics of the aerodynamic forces acting on the tower.In the following sections, the results for cases λ_3, λ_6, and λ_10 are discussed.

Pressure Coefficients on the Tower
Figure 7 shows the distributions of Cp on the tower surface at the mid-heights of each tower segment.When λ = 10, the distributions of Cp around the tower are almost symmetric with respect to θ = 0 for all segments.However, when λ = 3 and λ = 6, the angle θ at which Cp is maximized shifts in the −θ direction at upper segments, particularly at segments 1 and 2. This shifting is caused by the diversion of the flow approaching the tower by the rotating blades.In both figures, at all values of Ψ, the blade near the tower significantly diverts the flow approaching the tower, and the stagnation point on the tower shifts in the −θ direction.The degree of the diversion of the flow is more significant when λ = 3 than when λ = 6.Conversely, when λ = 10, the degree of the diversion of the flow due to the blade near the tower is very small at all values of Ψ, as shown in Figure 8c.
Energies 2017, 10, 121 9 of 18 flow approaching the tower, and the stagnation point on the tower shifts in the −θ direction.The degree of the diversion of the flow is more significant when λ = 3 than when λ = 6.Conversely, when λ = 10, the degree of the diversion of the flow due to the blade near the tower is very small at all values of Ψ, as shown in Figure 8c.Hence, the flow diversion caused by the rotating blades becomes less significant as λ increases, which could be due to the decrease of the aerodynamic force acting on the blades in the −Ψ direction arising from the reduction of the angle of attack on the blade (α, which is defined in arising from the reduction of the angle of attack on the blade (α, which is defined in Figure 1c), as shown in Table 2.This leads to smaller drag and lift coefficients of the blades.Previous studies [19,23] indicate that the lift coefficient of the S826 airfoil is 0 at α ≈ −5 • and is maximized at α ≈ 15 • .Additionally, the drag coefficient of the S826 airfoil is in the range of approximately 0.5-1.0 when −10 • ≤ α ≤ 15 • and increases sharply at α ≥ 15 • owing to the onset of a stall.The decreases in α and the aerodynamic force acting on the blades with an increase in λ leads to higher u values in the region between the rotor and tower, as shown in Figure 8a-c; higher p values on the downwind side of the blade, as shown in Figure 9; and higher maximum C p values around the tower at the mid-height of segment 1, as shown in Figure 7. Figure 1c), as shown in Table 2.This leads to smaller drag and lift coefficients of the blades.
Previous studies [19,23] indicate that the lift coefficient of the S826 airfoil is 0 at α ≈ −5° and is maximized at α ≈ 15°.Additionally, the drag coefficient of the S826 airfoil is in the range of approximately 0.5-1.0 when −10° ≤ α ≤ 15° and increases sharply at α ≥ 15° owing to the onset of a stall.The decreases in α and the aerodynamic force acting on the blades with an increase in λ leads to higher u values in the region between the rotor and tower, as shown in Figure 8a-c; higher p values on the downwind side of the blade, as shown in Figure 9; and higher maximum Cp values around the tower at the mid-height of segment 1, as shown in Figure 7.When λ = 3, the maximum values of Cp at the mid-heights of Segments 3-5 are larger than those of Segments 1 and 2, as shown in Figure 7a.This is mainly attributed to the decrease in the drag force acting on the blades, which reduces the wind velocity.The decrease in the drag force acting on the blade with an increase in the radial distance from the rotational axis (or with a decrease in the distance from the floor when the blade is at Ψ = 0°) is mainly attributed to the decrease in α, as shown in Table 2.When λ = 6 and λ = 10, the maximum values of Cp at the mid-heights of Segments 3-5 are smaller than those of Segments 1 and 2, as shown in Figure 7b,c, respectively.This tendency is more significant when λ = 10, mainly because of the increase in the lift force acting on the blade, which reduces the wind velocity.The increase in the lift force acting on the blades with an increase in the radial distance from the rotational axis could be caused by the increase in the relative wind velocity experienced by the blade.When λ = 10, the increase in the lift force acting on the blades with an increase in the radial distance from the rotational axis could also be a result of the increase in the lift coefficient of the blade with an increase in α.
In Figure 8d-f, the flow approaching the tower is significantly diverted by the blade when Ψ = 0°.However, for all other values of Ψ, the diversion of the flow approaching the tower by the blades is very small.This explains that the amount of shifting of θ at which Cp is maximized is very small at the mid-heights of Segments 3-5 for all the λ values shown in Figure 7.
Figure 10 shows the distributions of Cp' on the tower surface at the mid-heights of each tower segment.For all λ values, the distributions of Cp' around the tower are almost symmetrical with respect to θ = 0 at the segments near the floor.However, with the exception of these segments, the values of Cp' are generally higher in the +θ region than in the −θ region.This tendency could be caused by two factors when a blade passes the tower segments.The first factor is the low pressure region that forms on the downwind side of the blade.This region lowers the pressure on the tower.When λ = 3, the maximum values of C p at the mid-heights of Segments 3-5 are larger than those of Segments 1 and 2, as shown in Figure 7a.This is mainly attributed to the decrease in the drag force acting on the blades, which reduces the wind velocity.The decrease in the drag force acting on the blade with an increase in the radial distance from the rotational axis (or with a decrease in the distance from the floor when the blade is at Ψ = 0 • ) is mainly attributed to the decrease in α, as shown in Table 2.When λ = 6 and λ = 10, the maximum values of C p at the mid-heights of Segments 3-5 are smaller than those of Segments 1 and 2, as shown in Figure 7b,c, respectively.This tendency is more significant when λ = 10, mainly because of the increase in the lift force acting on the blade, which reduces the wind velocity.The increase in the lift force acting on the blades with an increase in the radial distance from the rotational axis could be caused by the increase in the relative wind velocity experienced by the blade.When λ = 10, the increase in the lift force acting on the blades with an increase in the radial distance from the rotational axis could also be a result of the increase in the lift coefficient of the blade with an increase in α.
In Figure 8d-f, the flow approaching the tower is significantly diverted by the blade when Ψ = 0 • .However, for all other values of Ψ, the diversion of the flow approaching the tower by the blades is very small.This explains that the amount of shifting of θ at which C p is maximized is very small at the mid-heights of Segments 3-5 for all the λ values shown in Figure 7.
Figure 10 shows the distributions of C p ' on the tower surface at the mid-heights of each tower segment.For all λ values, the distributions of C p ' around the tower are almost symmetrical with respect to θ = 0 at the segments near the floor.However, with the exception of these segments, the values of C p ' are generally higher in the +θ region than in the −θ region.This tendency could be caused by two factors when a blade passes the tower segments.The first factor is the low pressure region that forms on the downwind side of the blade.This region lowers the pressure on the tower.The second factor is the diversion of the flow approaching the tower by the blade, which shifts the stagnation point on the tower in the -θ direction.
Energies 2017, 10, 121 11 of 18 The second factor is the diversion of the flow approaching the tower by the blade, which shifts the stagnation point on the tower in the -θ direction.Figure 11 shows the pressure distribution at the mid-height of segment 4. As indicated by the figure, the pressure near the tower is reduced for all λ values as the blade approaches the tower.The pressure recovers when the blade moves away from the tower.Figure 12 shows the distributions of the instantaneous pressure coefficient Cpi (=2(p ) on the tower surface at the mid-heights of segment 4 when the blade nearest the tower is at various Ψ values.As shown in this figure, the θ at which Cpi is maximized on the upwind side of the tower generally shifts in theθ direction when the blade is near the tower (such as when Ψ = -20° to 20°) for all λ values.This shift in θ at which Cpi is maximized is due to the shifting of the stagnation point on the tower.When the blade approaches the tower, the values of Cpi on the +θ side are generally lower than those on the −θ-side.This is because the distance from the lower-pressure region on the downwind side of the blade is shorter on the +θ side than that on the -θ-side.Conversely, when the blade moves away from the tower, the values of Cpi on the -θ-side are generally lower than those on the +θ side.This is due to the distance from the lower-pressure region on the downwind side of the blade.However, when the blade moves away from the tower, the decrease in p on the −θ-side surface of the tower is suppressed by the shifting of θ at which Cpi is maximized in the −θ direction owing to the shifting of the stagnation point.Therefore, the degree of decrease in p on the −θ-side surface of the tower when the blade moves away from the tower is generally smaller than that on the +θ-side surface of the tower when the blade approaches the tower.As a result, the values of Cp' are generally higher on the +θ-side surface of the tower than those on the −θ-side surface.Figure 11 shows the pressure distribution at the mid-height of segment 4. As indicated by the figure, the pressure near the tower is reduced for all λ values as the blade approaches the tower.The pressure recovers when the blade moves away from the tower.Figure 12 shows the distributions of the instantaneous pressure coefficient C pi (=2 p − p ref /ρU ref 2 ) on the tower surface at the mid-heights of segment 4 when the blade nearest the tower is at various Ψ values.As shown in this figure, the θ at which C pi is maximized on the upwind side of the tower generally shifts in the -θ direction when the blade is near the tower (such as when Ψ = -20 • to 20 • ) for all λ values.This shift in θ at which C pi is maximized is due to the shifting of the stagnation point on the tower.When the blade approaches the tower, the values of C pi on the +θ side are generally lower than those on the −θ-side.This is because the distance from the lower-pressure region on the downwind side of the blade is shorter on the +θ side than that on the -θ-side.Conversely, when the blade moves away from the tower, the values of C pi on the -θ-side are generally lower than those on the +θ side.This is due to the distance from the lower-pressure region on the downwind side of the blade.However, when the blade moves away from the tower, the decrease in p on the −θ-side surface of the tower is suppressed by the shifting of θ at which C pi is maximized in the −θ direction owing to the shifting of the stagnation point.Therefore, the degree of decrease in p on the −θ-side surface of the tower when the blade moves away from the tower is generally smaller than that on the +θ-side surface of the tower when the blade approaches the tower.As a result, the values of C p ' are generally higher on the +θ-side surface of the tower than those on the −θ-side surface.

Drag and Lift Coefficients of the Tower
Figure 13 shows the vertical distributions of C Fx and C Fy of each tower segment.Figure 13a indicates that C Fx is maximized at Segment 6 when λ = 3 and λ = 6, and at Segment 7 when λ = 10, in Figure 13a.Hence, C Fx is the maximum at a height close to and lower than the lowest point of the rotor for all values of λ.This is because at Segment 6 or 7, the values of p are the lowest or second-lowest downwind of the tower, whereas the values of p upwind of the tower are approximately equal to those at the lower segments, as shown in Figure 7.

Drag and Lift Coefficients of the Tower
Figure 13 shows the vertical distributions of CFx and CFy of each tower segment.Figure 13a indicates that CFx is maximized at Segment 6 when λ = 3 and λ = 6, and at Segment 7 when λ = 10, in Figure 13a.Hence, CFx is the maximum at a height close to and lower than the lowest point of the rotor for all values of λ.This is because at Segment 6 or 7, the values of p are the lowest or second-lowest downwind of the tower, whereas the values of p upwind of the tower are approximately equal to those at the lower segments, as shown in Figure 7.The aerodynamic force acting on the upper segments of the tower could be very important from the perspective of the overturning moment of the tower.In Figure 13a, the values of CFx at Segments 1 and 2 increase as λ increases.This is because the values of p upwind of the tower increase as λ increases as shown in Figure 7.The values of CFy shown in Figure 13b are very small compared with those of CFx.In Figure 13b, the values of CFy at Segments 1 and 2 increase as λ decreases.This is mainly because the values of p on the -θ-side increase owing to the increased shift in the stagnation point in the −θ direction as λ decreases, as shown in Figure 7.The aerodynamic force acting on the upper segments of the tower could be very important from the perspective of the overturning moment of the tower.In Figure 13a, the values of C Fx at Segments 1 and 2 increase as λ increases.This is because the values of p upwind of the tower increase as λ increases as shown in Figure 7.The values of C Fy shown in Figure 13b are very small compared with those of C Fx .In Figure 13b, the values of C Fy at Segments 1 and 2 increase as λ decreases.This is mainly because the values of p on the -θ-side increase owing to the increased shift in the stagnation point in the −θ direction as λ decreases, as shown in Figure 7.

Drag and Lift Coefficients of the Tower
Figure 13 shows the vertical distributions of CFx and CFy of each tower segment.Figure 13a indicates that CFx is maximized at Segment 6 when λ = 3 and λ = 6, and at Segment 7 when λ = 10, in Figure 13a.Hence, CFx is the maximum at a height close to and lower than the lowest point of the rotor for all values of λ.This is because at Segment 6 or 7, the values of p are the lowest or second-lowest downwind of the tower, whereas the values of p upwind of the tower are approximately equal to those at the lower segments, as shown in Figure 7.The aerodynamic force acting on the upper segments of the tower could be very important from the perspective of the overturning moment of the tower.In Figure 13a, the values of CFx at Segments 1 and 2 increase as λ increases.This is because the values of p upwind of the tower increase as λ increases as shown in Figure 7.The values of CFy shown in Figure 13b are very small compared with those of CFx.In Figure 13b, the values of CFy at Segments 1 and 2 increase as λ decreases.This is mainly because the values of p on the -θ-side increase owing to the increased shift in the stagnation point in the −θ direction as λ decreases, as shown in Figure 7.  Figure 14 shows the vertical distributions of C Fx ' and C Fy ' of each tower segment.Here, the values of C Fx ' are smaller than those of C Fy ' at all segments for all values of λ. Figure 14a shows that for all values of λ, the values of C Fx ' at Segments 1-5, which are higher than the lowest point of the rotor, are larger than the values of C Fx ' at the segments near the floor.
Additionally, the values of C Fx ' are particularly large at Segments 3-5 and are maximized at one of these segments for all values of λ.The large fluctuations in F x at these segments are caused by the decrease and subsequent increase in F x that occur when a blade approaches and then moves away from the tower, as shown in Figure 15.This figure shows the time-history curves for the instantaneous drag coefficient C Fxi (=2F x /ρhdU ref 2 ) of each tower segment.
Energies 2017, 10, 121 14 of 18 Figure 14 shows the vertical distributions of CFx' and CFy' of each tower segment.Here, the values of CFx' are smaller than those of CFy' at all segments for all values of λ. Figure 14a shows that for all values of λ, the values of CFx' at Segments 1-5, which are higher than the lowest point of the rotor, are larger than the values of CFx' at the segments near the floor.
Additionally, the values of CFx' are particularly large at Segments 3-5 and are maximized at one of these segments for all values of λ.The large fluctuations in Fx at these segments are caused by the decrease and subsequent increase in Fx that occur when a blade approaches and then moves away from the tower, as shown in Figure 15.This figure shows the time-history curves for the instantaneous drag coefficient CFxi (=2Fx ρhdU ref 2

⁄
) of each tower segment.The decrease and subsequent increase in F x at Segments 3-5 that occur when a blade approaches and then moves away from the tower are mainly due to the decrease in the values of p upwind of the tower segments caused by the lower-pressure region formed downwind of the blade, as discussed in the previous section.With an increase in λ, the magnitude of the decrease in F x increases, as shown in Figure 15, and the maximum values of C Fx ' at Segments 3-5 increase, as shown in Figure 14a, because the values of p downwind of the blade decrease.The value of C Fx ' at Segment 1 in Figure 14a for λ = 10 is smaller than that for λ = 3 and λ = 6.This is because when λ = 10, the decrease in p downwind of the blade that passes the tower is very small owing to the small α.The values of C Fx ' at λ = 3 and λ = 6 are approximately equal.However, the fluctuations in p are larger upwind of the tower when λ = 3 and downwind of the tower when λ = 6, as shown in Figure 10.
Figure 14b shows that like C Fx ', the values of C Fy ' at Segments 1-5 are larger than the values of C Fy ' at the segments near the floor for all values of λ.When λ = 3 and λ = 6, the values of C Fy ' are maximized at Segment 4. When λ = 10, the value of C Fy ' at segment 4 is the second-largest.At Segment 6, which is lower than the lowest point of the rotor, the value of C Fy ' is maximized.The values of C Fy ' are relatively large at Segments 3-5 for all values of λ because when a blade approaches the tower, F y increases and then decreases, and when the blade moves away from the tower, F y decreases further and then increases, as shown in Figure 16. Figure 16 shows the time history of the instantaneous lift coefficient C Fyi (=2F y /ρhdU ref 2 ) of each tower segment.Periodic fluctuation in F y occurs because the low-pressure region downwind of the blade has a greater effect on the pressure on the +θ-side surface of the tower segments when the blade approaches the tower and a greater effect on the pressure on the −θ-side surface of the tower segments when the blade moves away from the tower, as shown in Figure 11.C Fy ' is maximized at Segment 6 when λ = 10 because in addition to the aforementioned fluctuation in F y that occurs when a blade passes the tower, the F y values fluctuate with the largest amplitude and smaller frequency than the blade-passing frequency, as shown in Figure 16.In this figure, F y is negative at t ≈ 2.72-2.77s and positive at t ≈ 2.77-2.81s, with a large amplitude.This is because a strong vortex tends to form on one lateral side of segment 6 for longer periods than the blade-passing interval, followed by a similar formation on the other lateral side.Figure 14b shows that at segment 1, C Fy ' is larger for λ = 6 than the values obtained when λ = 3 and λ = 10.This is because the value of p on the tower segment in the vicinity of θ = ± 90 • fluctuates more when λ = 6, as shown in Figure 10.
Energies 2017, 10, 121 15 of 18 The decrease and subsequent increase in Fx at Segments 3-5 that occur when a blade approaches and then moves away from the tower are mainly due to the decrease in the values of p upwind of the tower segments caused by the lower-pressure region formed downwind of the blade, as discussed in the previous section.With an increase in λ, the magnitude of the decrease in Fx increases, as shown in Figure 15, and the maximum values of CFx' at Segments 3-5 increase, as shown in Figure 14a, because the values of p downwind of the blade decrease.The value of CFx' at Segment 1 in Figure 14a for λ = 10 is smaller than that for λ = 3 and λ = 6.This is because when λ = 10, the decrease in p downwind of the blade that passes the tower is very small owing to the small α.The values of CFx' at λ = 3 and λ = 6 are approximately equal.However, the fluctuations in p are larger upwind of the tower when λ = 3 and downwind of the tower when λ = 6, as shown in Figure 10.
Figure 14b shows that like CFx', the values of CFy' at Segments 1-5 are larger than the values of CFy' at the segments near the floor for all values of λ.When λ = 3 and λ = 6, the values of CFy' are maximized at Segment 4. When λ = 10, the value of CFy' at segment 4 is the second-largest.At Segment 6, which is lower than the lowest point of the rotor, the value of CFy' is maximized.The values of CFy' are relatively large at Segments 3-5 for all values of λ because when a blade approaches the tower, Fy increases and then decreases, and when the blade moves away from the tower, Fy decreases further and then increases, as shown in Figure 16. Figure 16 shows the time history of the instantaneous lift coefficient CFyi (=2Fy ρhdU ref 2

⁄
) of each tower segment.Periodic fluctuation in Fy occurs because the low-pressure region downwind of the blade has a greater effect on the pressure on the +θ-side surface of the tower segments when the blade approaches the tower and a greater effect on the pressure on the −θ-side surface of the tower segments when the blade moves away from the tower, as shown in Figure 11.
CFy' is maximized at Segment 6 when λ = 10 because in addition to the aforementioned fluctuation in Fy that occurs when a blade passes the tower, the Fy values fluctuate with the largest amplitude and smaller frequency than the blade-passing frequency, as shown in Figure 16.In this figure, Fy is negative at t ≈ 2.72-2.77s and positive at t ≈ 2.77-2.81s, with a large amplitude.This is because a strong vortex tends to form on one lateral side of segment 6 for longer periods than the blade-passing interval, followed by a similar formation on the other lateral side.Figure 14b shows that at segment 1, CFy' is larger for λ = 6 than the values obtained when λ = 3 and λ = 10.This is because the value of p on the tower segment in the vicinity of θ = ± 90° fluctuates more when λ = 6, as shown in Figure 10.

Conclusions
As a first step toward clarifying the effects of aerodynamic forces acting on the tower of a wind turbine on the vibratory characteristics of the tower, we investigated the effects of the rotating blades of an upwind-type three-blade HAWT on the basic characteristics of the aerodynamic forces acting on its tower.This was done by conducting improved delayed detached-eddy simulations of the wind flow around the tower and the wind turbine operating at tip-speed ratios of λ = 3 (low), λ = 6 (optimum), and λ = 10 (high).The tower was divided into 10 equal segments, and the aerodynamic forces acting on each segment were analyzed.
We discovered that both the diversion of the flow approaching the tower by the rotating blades and the low-pressure region that forms downwind of the blades have significant effects on the aerodynamic forces acting on the tower.For example, owing to the diversion of the flow by the blades, the azimuth angle around the tower at which the pressure reaches a maximum at each height can shift significantly in the direction of the movement of the blade passing the tower.In addition, because of the low-pressure region downwind of the blades, the fluctuation in the lift force of the tower can be significantly larger than that in its drag force.
These characteristics of the aerodynamic forces acting on the tower might contribute significantly to the occurrence of unacceptably large stresses and deflections in the side-to-side direction when the blade passing frequency coincides with the natural frequency of the tower.Owing to the negligible contribution of the aerodynamic damping of HAWTs in the side-to-side direction, the damping ratio for the side-to-side oscillations of the tower is generally one order of magnitude less than that for the fore-aft oscillations of the tower [2].Therefore, in future research, we plan to investigate the effects of the aerodynamic forces acting on the tower that are influenced by the rotating blades on the vibratory characteristics of the tower in the side-to-side direction by

Conclusions
As a first step toward clarifying the effects of aerodynamic forces acting on the tower of a wind turbine on the vibratory characteristics of the tower, we investigated the effects of the rotating blades of an upwind-type three-blade HAWT on the basic characteristics of the aerodynamic forces acting on its tower.This was done by conducting improved delayed detached-eddy simulations of the wind flow around the tower and the wind turbine operating at tip-speed ratios of λ = 3 (low), λ = 6 (optimum), and λ = 10 (high).The tower was divided into 10 equal segments, and the aerodynamic forces acting on each segment were analyzed.
We discovered that both the diversion of the flow approaching the tower by the rotating blades and the low-pressure region that forms downwind of the blades have significant effects on the aerodynamic forces acting on the tower.For example, owing to the diversion of the flow by the blades, the azimuth angle around the tower at which the pressure reaches a maximum at each height can shift significantly in the direction of the movement of the blade passing the tower.In addition, because of the low-pressure region downwind of the blades, the fluctuation in the lift force of the tower can be significantly larger than that in its drag force.
These characteristics of the aerodynamic forces acting on the tower might contribute significantly to the occurrence of unacceptably large stresses and deflections in the side-to-side direction when the blade passing frequency coincides with the natural frequency of the tower.Owing to the negligible contribution of the aerodynamic damping of HAWTs in the side-to-side direction, the damping ratio for the side-to-side oscillations of the tower is generally one order of magnitude less than that for the fore-aft oscillations of the tower [2].Therefore, in future research, we plan to investigate the effects of

Figure 1 .
Figure 1.Geometry of the wind turbine and tower and the definitions of the angles Ψ, α, γ, and θ.(a) Side view of the wind turbine and tower; (b) front view of the wind turbine and tower; (c) cross-sectional view of the blade at Ψ = 0°; and (d) horizontal cross-sectional view of the tower.

Figure 1 .
Figure 1.Geometry of the wind turbine and tower and the definitions of the angles Ψ, α, γ, and θ.(a) Side view of the wind turbine and tower; (b) front view of the wind turbine and tower; (c) cross-sectional view of the blade at Ψ = 0 • ; and (d) horizontal cross-sectional view of the tower.

Figure 2 .
Figure 2. Computational grid and boundary conditions.(a) Entire computational domain; (b) rotational domain; and (c) horizontal cross-section of the grid near the tower at the mid-height of Segment 3.

Figure 2 .
Figure 2. Computational grid and boundary conditions.(a) Entire computational domain; (b) rotational domain; and (c) horizontal cross-section of the grid near the tower at the mid-height of Segment 3.

Figure 3
Figure 3 compares the values of the power coefficient C P (=2QΩ/ρAU ref 3 ) and the thrust coefficient C T (=2T/ρAU ref 2 ) evaluated from the simulations to those obtained from the wind-tunnel

Figure 3 .
Figure 3. Computational and experimental results for the wind-turbine performance.(a) Power coefficient; and (b) thrust coefficient.

Figure 4
Figure4compares the simulation and wind-tunnel experiment results for the distribution in the y direction of at the hub height in the wind-turbine wake at x/D = 1.The simulated values for strongly match the experimental results[17,18] qualitatively and quantitatively.The discrepancies of between case λ_6 and other cases with λ = 6 are very small.

Figure 5 2 ⁄
Figure 5 shows the sensitivity of DR, the grid resolution of the structured mesh region near the tower, and dt on the vertical distributions of the mean drag coefficient CFx (=2F x ρhdU ref 2 ), the mean lift coefficient CFy (=2F y ρhdU ref 2 ), the fluctuating drag coefficient CFx'(=2σFx ρhdU ref 2 ⁄ ) and the fluctuating lift coefficient CFy'(=2σFy ρhdU ref 2 ⁄) of each tower segment.Here, Fx and Fy are the drag and lift forces, respectively, acting on each tower segment, and σFx and σFy denote the root mean squares of Fx and Fy, respectively.Except for case G_coarse, which is not shown in the figures, the discrepancies of these drag and lift coefficients between case λ_6 and other cases with λ = 6 are small.

Figure 3 .
Figure 3. Computational and experimental results for the wind-turbine performance.(a) Power coefficient; and (b) thrust coefficient.

Figure 4
Figure4compares the simulation and wind-tunnel experiment results for the distribution in the y direction of u at the hub height in the wind-turbine wake at x/D = 1.The simulated values for u strongly match the experimental results[17,18] qualitatively and quantitatively.The discrepancies of u between case λ_6 and other cases with λ = 6 are very small.

Figure 3 .
Figure 3. Computational and experimental results for the wind-turbine performance.(a) Power coefficient; and (b) thrust coefficient.

Figure 4
Figure4compares the simulation and wind-tunnel experiment results for the distribution in the y direction of at the hub height in the wind-turbine wake at x/D = 1.The simulated values for strongly match the experimental results[17,18] qualitatively and quantitatively.The discrepancies of between case λ_6 and other cases with λ = 6 are very small.

Figure 5
Figure 5 shows the sensitivity of DR, the grid resolution of the structured mesh region near the tower, and dt on the vertical distributions of the mean drag coefficient CFx (=2F x ρhdU ref 2 ), the mean lift coefficient CFy (=2F y ρhdU ref 2 ), the fluctuating drag coefficient CFx'(=2σFx ρhdU ref 2 ⁄ ) and the fluctuating lift coefficient CFy'(=2σFy ρhdU ref 2 ⁄) of each tower segment.Here, Fx and Fy are the drag and lift forces, respectively, acting on each tower segment, and σFx and σFy denote the root mean squares of Fx and Fy, respectively.Except for case G_coarse, which is not shown in the figures, the discrepancies of these drag and lift coefficients between case λ_6 and other cases with λ = 6 are small.

Figure 5 18 Figure 5 .
Figure 5 shows the sensitivity of D R , the grid resolution of the structured mesh region near the tower, and dt on the vertical distributions of the mean drag coefficient C Fx (=2F x /ρhdU ref 2 ), the mean lift coefficient C Fy (=2F y /ρhdU ref 2 ), the fluctuating drag coefficient C Fx '(=2σ Fx /ρhdU ref 2 ) and the fluctuating lift coefficient C Fy '(=2σ Fy /ρhdU ref 2 ) of each tower segment.Here, F x and F y are the drag and lift forces, respectively, acting on each tower segment, and σ Fx and σ Fy denote the root mean

Figure 6 .
Figure 6.Sensitivity of DR, the grid resolution around the tower, and dt on distributions of the mean and fluctuating pressure coefficients on the tower surface at the mid-height of Segment 6 at λ = 6.(a) Mean pressure coefficient; and (b) fluctuating pressure coefficient.

Figure 5 .
Figure 5. Sensitivity of D R , the grid resolution near tower, and dt on vertical distributions of the mean and fluctuating drag and lift coefficients of each tower segment at λ = 6.(a) Mean drag coefficient; (b) mean lift coefficient; (c) fluctuating drag coefficient; and (d) fluctuating lift coefficient.

Figure 5 .
Figure 5. Sensitivity of DR, the grid resolution near tower, and dt on vertical distributions of the mean and fluctuating drag and lift coefficients of each tower segment at λ = 6.(a) Mean drag coefficient; (b) mean lift coefficient; (c) fluctuating drag coefficient; and (d) fluctuating lift coefficient.

Figure 6 .
Figure 6.Sensitivity of DR, the grid resolution around the tower, and dt on distributions of the mean and fluctuating pressure coefficients on the tower surface at the mid-height of Segment 6 at λ = 6.(a) Mean pressure coefficient; and (b) fluctuating pressure coefficient.

Figure 6 .
Figure 6.Sensitivity of D R , the grid resolution around the tower, and dt on distributions of the mean and fluctuating pressure coefficients on the tower surface at the mid-height of Segment 6 at λ = 6.(a) Mean pressure coefficient; and (b) fluctuating pressure coefficient.

Energies 2017, 10 , 121 8 of 18 Figure 6
Figure 6 shows the sensitivity of DR, the grid resolution of the structured mesh region near the tower, and dt on the horizontal distributions of the mean pressure coefficient Cp (=2(p − p ref ) ρU ref 2 ⁄ )

FigureFigure 7 .
Figure 8a,b show the streamlines at the mid-height of Segment 1 for λ = 3 and λ = 6, respectively.In both figures, at all values of Ψ, the blade near the tower significantly diverts the Figure 7. Distributions of the mean pressure coefficient on the tower surface at the mid-heights of each tower segment.(a) λ = 3; (b) λ = 6; and (c) λ = 10.

Figure
Figure 8a,b show the streamlines at the mid-height of Segment 1 for λ = 3 and λ = 6, respectively.In both figures, at all values of Ψ, the blade near the tower significantly diverts the flow approaching the tower, and the stagnation point on the tower shifts in the −θ direction.The degree of the diversion of the flow is more significant when λ = 3 than when λ = 6.Conversely, when λ = 10, the degree of the diversion of the flow due to the blade near the tower is very small at all values of Ψ, as shown in Figure8c.

Figure 12 .
Figure 12.Distributions of the instantaneous pressure coefficient on the tower surface at the mid-height of Segment 4 for moments when the blade nearest the tower is at various Ψ values.(a) λ = 3; (b) λ = 6; and (c) λ = 10.

Figure 12 .
Figure 12.Distributions of the instantaneous pressure coefficient on the tower surface at the mid-height of Segment 4 for moments when the blade nearest the tower is at various Ψ values.(a) λ = 3; (b) λ = 6; and (c) λ = 10.

Figure 12 .
Figure 12.Distributions of the instantaneous pressure coefficient on the tower surface at the mid-height of Segment 4 for moments when the blade nearest the tower is at various Ψ values.(a) λ = 3; (b) λ = 6; and (c) λ = 10.

Figure 13 .
Figure 13.Vertical distributions of the mean drag and mean lift coefficients of each tower segment.(a) Mean drag coefficient; and (b) mean lift coefficient.

Figure 14 .
Figure 14.Vertical distributions of the fluctuating drag and fluctuating lift coefficients of each tower segment.(a) Fluctuating drag coefficient; and (b) fluctuating lift coefficient.

Figure 13 .
Figure 13.Vertical distributions of the mean drag and mean lift coefficients of each tower segment.(a) Mean drag coefficient; and (b) mean lift coefficient.

Figure 13 .
Figure 13.Vertical distributions of the mean drag and mean lift coefficients of each tower segment.(a) Mean drag coefficient; and (b) mean lift coefficient.

Figure 14 .
Figure 14.Vertical distributions of the fluctuating drag and fluctuating lift coefficients of each tower segment.(a) Fluctuating drag coefficient; and (b) fluctuating lift coefficient.

Figure 14 .
Figure 14.Vertical distributions of the fluctuating drag and fluctuating lift coefficients of each tower segment.(a) Fluctuating drag coefficient; and (b) fluctuating lift coefficient.

Figure 15 .
Figure 15.Time history of the instantaneous drag coefficient of each tower segment.The vertical dotted lines indicate the time at which the blade nearest the tower is at Ψ = 0°.(a) λ = 3; (b) λ = 6; and (c) λ = 10.

Figure 15 .
Figure 15.Time history of the instantaneous drag coefficient of each tower segment.The vertical dotted lines indicate the time at which the blade nearest the tower is at Ψ = 0 • .(a) λ = 3; (b) λ = 6; and (c) λ = 10.

Figure 16 .
Figure 16.Time history of the instantaneous lift coefficient of each tower segment.The vertical dotted lines indicate the time at which the blade nearest the tower is at Ψ = 0°.(a) λ = 3; (b) λ = 6; and (c) λ = 10.

Figure 16 .
Figure 16.Time history of the instantaneous lift coefficient of each tower segment.The vertical dotted lines indicate the time at which the blade nearest the tower is at Ψ = 0 • .(a) λ = 3; (b) λ = 6; and (c) λ = 10.

Table 1 .
Computational conditions.Here, λ is the tip-speed ratio; D R is the diameter of the rotational domain; D is the diameter of the rotor; d tr is the shortest distance between the tip of a blade and the edge of the rotational domain; and dt is the time-step size.

Table 2 .
Angle of attack of the blade at Ψ = 0 • at the mid-height of each tower segment.

Table 2 .
Angle of attack of the blade at Ψ = 0° at the mid-height of each tower segment.