**1. Introduction**

The adoption and promotion of renewable energy has gained widespread interest worldwide, including in Japan. For this reason, the introduction of wind farms with wind turbines has accelerated. However, there have also been grea<sup>t</sup> concerns. Unfortunately, an increasing trend in the number of serious accidents has been reported, such as the accidental fall of a wind turbine nacelle, particularly in wind power stations built on complex terrain in mountainous areas. Utmost caution is required, especially for wind power stations with complicated topography, in planning an optimum arrangemen<sup>t</sup> of wind turbines and controlling both maintenance and management. Diversified tasks are increasingly important to avoid accumulated fatigue of wind load in wind turbines due to terrain-induced turbulence, to reduce malfunctions and accidents inside and outside wind turbines, and to improve the availability of wind turbines [1–11]. Considering social and engineering requests at present, the crucial purpose of the present study is to establish a system with a numerical diagnostic technique for wind status, which contributes to the proper operation of wind farms, an adequate understanding of the indigenous wind environments of each site, including terrain-induced turbulence, and a reduction of malfunctions and accidents associated with wind turbines [12–21].

We conducted research with a demonstration not only to examine the impacts of terrain-induced turbulence on flapwise fatigue damage of wind turbine blades, but also to clarify the mechanism of terrain-induced turbulence generation. Specifically, we targeted a large-scale wind turbine built on complex mountainous terrain. Electric strain gauges were installed at the base of three blades of the wind turbine. We developed a measurement system to automatically obtain time-series data of strain fluctuation via the gauges. Subsequently, the damage equivalent load (DEL) [22,23] of the flapwise bending and vibration of wind turbine blades was calculated on the basis of the collected data. Furthermore, we examined the relationship with the output results of wind direction and speed sensors, which were installed on the wind turbine nacelle for the purpose of wind turbine control.

Firstly, the time period of the maximum DEL on the blades during the measurement period in this study was identified, as well as the time period of the maximum load regarding fatigue damage of the blades, and the situation was quantitatively examined. Simultaneously, a high-resolution simulation of numerical wind conditions was performed on the basis of large-eddy simulation (LES) to closely examine three-dimensional airflow structures when the wind turbine blades experienced terrain-induced turbulence. Secondly, the present study derived a formula, a relational expression indicating a correlation between the DEL calculated on the basis of actual measured data and the values of nacelle wind speed, which were obtained by the nacelle anemometer, that is, the fatigue damage of the blades. Using this formula and measurement data of wind conditions over a one-year period, the influence of terrain-induced turbulence on the age-related degradation of wind turbine blades was quantitatively assessed.

Lastly, through the analysis of the series of measurement data and numerical wind condition simulations, we proposed two types of new assessment scales to evaluate the impacts of terrain-induced turbulence on wind turbine blades and demonstrated the operation methods of these generalization indexes.

#### **2. Overview of the Kushikino Reimei Wind Farm**

The present study was conducted, with cooperation of Kyudenko New Energy Co., Ltd.(Fukuoka, Japan), focusing on the Kushikino Reimei Wind Farm (established in November 2012) in Hashima, Ichikikushikino City, Kagoshima Prefecture, Japan (Figure 1). This wind farm was installed with ten 2 MW (megawatt) downwind wind turbines (Hitachi, Ltd. (Tokyo, Japan)). The target wind turbine was wind turbine #10. The present study focused on the impacts of terrain-induced turbulence on the blades of wind turbine #10. Turbulence was generated when easterly winds passed over Mt. Benzaiten (elevation 519 m) (see Figures 2 and 3, and Table 1). Figure 4 provides an outline of the wind turbine. Figure 5 illustrates the power curve of the wind turbine. Figure 6 shows a vane anemometer which was installed on the wind turbine nacelle. The pitch and yaw control of the wind turbine was performed based on the sensor information, as shown in Figure 5. The present study analyzed the airflow field generated around wind turbine #10, utilizing output results of wind direction and speed sensors (research data results are described later).

**Figure 1.** Map of the Kushikino Reimei Wind Farm and the surrounding area.

**Figure 2.** Photo of wind turbine #10.

**Figure 3.** Relative locations of Mt. Benzaiten (elevation 519 m) and wind turbine #10.

**Table 1.** Elevation information for wind turbine #10 and distance between Mt. Benzaiten (elevation 519 m) and wind turbine #10.

**Figure 4.** Outline of a wind turbine.

**Figure 5.** Power curve of a wind turbine.

**Figure 6.** Nacelle propeller-vane anemometer.

#### **3. In-Situ Data Analysis**

#### *3.1. Analysis of Wind Turbine Power Output Data*

The theoretical power output was compared with the measured data for wind turbine #10. The results are shown in Figure 7. Figure 7a,b show measured values of northerly and easterly winds, respectively; the number of measured data values was 1578 for the northerly wind, obtained in 10-min periods in January of 2013, and similarly, 601 data values were obtained for the easterly wind in June of 2013. The easterly wind pattern shows a high dispersion compared with that of the northerly wind.

**Figure 7.** Comparison between theoretical data and measured data for wind turbine #10. (**a**) Northerly wind. (**b**) Easterly wind.

#### *3.2. Analysis of Wind Turbine Alarm Data*

A variety of alarm conditions exist for a wind turbine. For this research, the following two items regarding wind conditions were our focus, and the operation status of wind turbine #10 was inspected. The target period of the data analysis was 14 months (November 2012–January 2014).

(1) Shutdown due to excessive yaw error (definition: malfunctions occur due to deviation of the nacelle direction and the anemoscope direction).

(2) Discordance in wind directions of sensors for wind direction and speed (definition: malfunctions occur due to disagreement in the values of the two sensors).

Table 2 provides the number of alarm occurrences due to the above alarm items for all wind directions. The mean values of other wind turbines are shown for a comparative examination. This table indicates that the alarm of wind turbine #10 occurred frequently due to the two items mentioned above. In Table 3, alarm occurrences of Table 2 are shown depending on the wind directions. Figure 8 shows a graph generated from Table 3. Table 3 and Figure 8 clearly indicate that wind turbine #10 had an extremely high number of alarm occurrences when wind blew from the east compared with other wind directions. This phenomenon suggested that wind turbine #10 was affected over time by unsteady wind direction fluctuations when easterly winds occurred.

**Table 2.** Number of alarm occurrences for wind conditions under all wind directions.


**Table 3.** Number of alarm occurrences due to wind conditions for each wind direction.


**Figure 8.** Result of graphing the values in Table 3. (**a**) Number of alarm occurrences for yaw misalignment. (**b**) Number of alarm occurrences for wind direction mismatch of the wind vane. Note: Wind turbine systems resort to shutting down the wind turbine if the yaw misalignment exceeds a threshold.

#### *3.3. Analysis of Nacelle Propeller-Vane Anemometer Data*

Wind energy is input via wind turbine blades. Thus, it is absolutely essential to monitor the behavior of the flapwise bending and vibration at the base of the blades for the assessment of the durability of the entire wind turbine. In the present study, two types of electric strain gauges were installed at the base of three blades on wind turbine #10 (the base: the position at a distance of approximately 1.3 m from the hub-connected point, referred to in Figure 4). A measurement synchronization system was developed to measure both of these values and basic information on the wind turbine operation (eight items: wind direction, wind speed and azimuth of the nacelle, pitch angle, rotational speed of the generator, active power of the power conditioning system (PCS) system, azimuth angle, and longitudinal acceleration of the nacelle). Actual measurement data were collected at 50 Hz (0.02 interval, 50 cycles per second) via this measurement synchronization system. The average values of the two sets of wind direction and speed sensors for wind turbine control, which were installed on the wind turbine nacelle shown in Figure 6, were used as measurement data regarding the wind direction and speed of the nacelle. The measurement period was from 3 November 2015, 0:00 a.m. JST to 17 March 2016, 7:00 a.m. JST.

Wind roses for all of the wind velocity classes are shown in Figure 9, where the above-measured data of wind direction and wind speed of the nacelle were classified into 16 wind directions. Further, the numerical data are shown in Table 4. Regarding frequency distribution, the highest frequency, 22.5%, was for the northerly wind (Total: 19237, N: 4331) and the frequency of the easterly wind was 4.4% (Total: 19237, E: 856). The mean speed of the easterly wind was lower than that of the northerly wind; the northerly wind mean speed was 6.1 m/s, while the easterly mean speed was 4.5 m/s.

**Figure 9.** Frequency distribution of the direction of the 10-min average wind (%): (**a**) wind rose (frequency distribution) and the average of the 10-min average wind speed observed for 16 wind directions (m/s): (**b**) wind speed by direction (wind measurement height: hub height (60 m), analysis period: 3 November 2015, 0:00 a.m. JST–17 March 2016, 7:00 a.m. JST).

**Table 4.** Frequency distribution of the direction of the 10-min average wind (%) and the average of the 10-min average wind speed observed for 16 directions (wind measurement height: hub height (60 m), analysis period: 3 November 2015, 0:00 a.m. JST–17 March 2016, 7:00 a.m. JST).


Data groups were analyzed for nacelle wind direction and nacelle wind speed, and the standard deviation was calculated for nacelle wind speed values on the basis of 10-min intervals for the data measurement period of wind turbine #10 from 3 November 2015, 0:00 a.m. JST to 17 March 2016, 7:00 a.m. JST. Furthermore, measurement data were classified into 12 wind directions and were arranged with analysis results of the damage equivalent load (DEL), which is described later (see Table 5). Figures 10 and 11 illustrate the analysis results of the wind velocity standard deviation and turbulence intensity for both northerly and easterly winds. To examine the status of wind turbine generation, analysis data targets corresponded to winds of 4 m/s or higher.


**Table 5.** Wind direction range and total number of data values.

Note: includes only data from 10-min periods with an average wind speed of 4 m/s (cut-in wind speed) or higher.

**Figure 10.** Relationship between the standard deviation and the average of the wind speed in 10-min periods for two wind directions (wind measurement height: hub-height (60 m), analysis period: 3 November 2015, 0:00 a.m. JST–17 March 2016, 7:00 a.m. JST). (**a**) Northerly wind. (**b**) Easterly wind.

**Figure 11.** Relationship between turbulence intensity and the average of the wind speed in 10-min periods for two wind directions (wind measurement height: hub-height (60 m), analysis period: 3 November 2015, 0:00 a.m. JST–17 March 2016, 7:00 a.m. JST). (**a**) Northerly wind. (**b**) Easterly wind.

Wind velocity standard deviations and turbulence intensity were calculated by formula (1):

$$TurbulenceIntensity(TI) = \frac{\sigma\_u}{\overline{u}} = \frac{\sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left(u\_i'\right)^2}}{\overline{u}} \tag{1}$$

where:

$$
\mu' = \mu(t) - \overline{\mu} \tag{2}
$$

In Figures 10 and 11, values of standard deviations of easterly wind velocity and values of their corresponding turbulence intensity were notably larger in comparison to those of the northerly wind. In addition, it was confirmed that the turbulence intensity of the easterly wind, in a wind velocity class 10 m/s and lower, frequently exceeded the turbulence intensity category of the International Electrotechnical Commission (IEC). It was inferred that the influence of Mt. Benzaiten (elevation 519 m) at a distance of 300 m from wind turbine #10 was a major cause of the increased turbulence intensity when easterly winds occurred. This is examined later, employing numerical wind condition simulations.

#### *3.4. Analysis of Wind Turbine Blade Strain Data*

As mentioned above, wind energy is input via wind turbine blades. Thus, it is absolutely essential to monitor the behaviors of flapwise bending and vibration at the base of blades for the assessment of the durability of the entire wind turbine. In the present study, electric strain gauges were installed at the base of three blades on wind turbine #10, as shown in Figure 4. Figure 12 presents the blade strain data of two types of winds in the measurement period. They were compared between two time-history data waveforms when the mean values of the strain gauge of nacelle wind speeds were approximately 9 m/s. The northerly wind showed the highest occurrence frequency, and the easterly wind showed notably large values in turbulence intensity, as was described previously. Data that are framed by a continuous line represent the time zones of two wind patterns, which showed approximately 9 m/s mean values of the nacelle wind velocity (Figure 12). A comparison of the results of measurement data between the easterly wind in Figure 12b and the northerly wind in Figure 12a reveals that variable amplitudes of strain gauges were notably large for the easterly wind, and thus, the blades of the target wind turbine #10 vibrated due to the large flapwise wind loads.

This study conducted a further quantitative examination of the strain data of wind turbine blades. From strain data, time history data of flapwise bending moments were extracted with the cooperation of a wind turbine manufacturer. Upon applying the rainflow counting algorithm [22,23], the damage equivalent load (DEL) was calculated (see Formula 3). DEL is the most commonly applied index in the wind power industry in arguments regarding the fatigue damage of wind turbines. In the present research, the obtained values of DEL were normalized by a design value for a wind velocity of 12 m/s, employing the aero-elastic analysis software Bladed.

$$DEL = \left\{ \frac{\sum\_{i=1}^{n} \left( F\_i^m \cdot n\_i \right)}{N} \right\}^{\frac{1}{m}} \tag{3}$$

where

*Fi* is the load of the *i-*th class of the fatigue load spectrum;

*ni* is the number of cycles in the *i-*th class of the fatigue load spectrum;

*N* is the equivalent of cycles;

*m* is the S-N (stress-number of cycles to failure) curve slope for relevant material.

*N* = 600 and *m* = 10 with fiber-reinforced plastic (FRP) blades in this study.

**Figure 12.** Blade strain data (blade flapwise bending raw data). (**a**) Northerly wind. (**b**) Easterly wind. Note: interval: 0.02 seconds, average wind speed: approx. 9 m/s.

Figure 13 shows time variations at 10-min intervals for nacelle wind velocity, its standard deviation and DEL. Similar to Figure 12, time zones where the 10-min mean value of nacelle wind velocity was approximately 9 m/s are the focus of this discussion. In the time period of the northerly wind (November 9, 2015, 19:30–19:40 p.m. JST), the wind velocity value was 9.4 m/s. Corresponding to this time period, the standard deviation was 1.3 m/s and DEL was 0.99. Contrarily, the time period for which 9 m/s was intended to be a mean speed was from November 13, 2015, 9:40 to 9:50 a.m. JST for the easterly wind, showing approximately 9.1 m/s as the 10-min mean values of the nacelle wind velocity. The value of standard deviation for the easterly wind is worth mentioning in particular. The mean value of the standard deviation of the nacelle wind velocity (November 13, 2015, 9:40–9:50 a.m. JST) was 2.3 m/s, which was approximately 1.8 times higher than that of the northerly wind. Accordingly, the DEL was also notable at 2.03, which was approximately twice as large as that of the northerly wind. Conclusively, the easterly wind status, which had the issue of terrain-induced turbulence, was proven to have distinctive differences in comparison with the northerly wind status, which was scarcely affected by terrain-induced turbulence, while it had the highest frequency of occurrence. Based on these findings, it was revealed that for easterly wind patterns, the standard deviation was large. Accordingly, the DEL was large, and thus terrain-induced turbulence directly influenced the turbine blades. The DEL value of 2.03 for the easterly wind (November 13, 2015, 9:40–9:50 a.m. JST) was the highest in the entire present study; that is, this period was the time period where the largest load was generated regarding the fatigue damage of the blades of the wind turbine. This scale of 2.03 in DEL means that if airflow with a property of DEL = 2.03 continued for 5.88 years, the total load on the wind turbine blades would reach the design load for the designed service life.

**Figure 13.** Time series of wind speed, standard deviation and normalized damage equivalent load (DEL) for flapwise blade bending. Plotted values were evaluated from instantaneous data values within 10-min periods. (**a**) Northerly wind. (**b**) Easterly wind.

#### **4. Numerical Simulation of Airflow with the WRF Mesoscale Model**

As described in Section 3.4, the wind vanes on the nacelle of WT#10 indicated the presence of an easterly wind from 9:40 to 9:50 a.m. JST, November 13, 2015. To investigate the wind flow from this time period closely, a numerical simulation was performed using the Weather Research and Forecasting (WRF) mesoscale model [24] in the present study. This mesoscale model was developed cooperatively by a number of institutions, including NCAR (National Center for Atmospheric Research), the University of Oklahoma, NCEP (US National Center for Environment Prediction), NOAA (Forecast System Laboratory of the National Oceanic and Atmospheric Administration), and AFWA (Air Force Weather Agency). WRF is a three-dimensional, fully compressible, non-hydrostatic model that has been used both operationally and for research worldwide. WRF is considered a successor to the non-hydrostatic MM5 model, which was developed principally by NCAR. A number of physics models are included in WRF, such as radiation models to calculate solar and atmospheric radiation, a turbulence model that simulates turbulence in the mixed-layer, a cloud physics model that takes into

account water vapor, cloud water, rainwater, snow, and hail, and a land surface model that calculates the surface temperature, soil temperature, soil water content, snowfall, and surface flux. Furthermore, WRF allows the use of the latest physics models and a data assimilation system; thus, WRF is suitable for predicting and simulating localized heavy rainfall, gusts, and other meteorological phenomena.

Figure 14 illustrates the computational domain used for the present simulation. In the present study, four-layer nesting was adopted. The spatial resolution of the smallest domain (d04), in which the wind farm is located, was 333.33 m. The height of the smallest grid cell was approximately 8 m. For the terrain data set, the model used the global digital elevation model GTOPO30 with spatial resolutions of 30 seconds in latitude and longitude provided by the United States Geological Survey (USGS) and higher-resolution data; specifically, the 50-m digital elevation model (GSI50, spatial resolution: 1.5 seconds in the latitude direction and 2.25 seconds in the longitude direction) provided by the Geospatial Information Authority of Japan (GSI). For land use and vegetation data, the USGS data that are pre-installed in WRF (spatial resolution: 1◦) were used. The meteorological GPV (grid point value) data that were used for the boundary conditions in the present study were the NCEP Final Analysis (NCEP-FNL) data (spatial resolution: 0.5◦; temporal resolution: six hours), which are global analysis data. These data were used every six hours (with no nudging). For sea surface temperature (SST) data, the skin temperature from the NCEP–FNL data was used. The airflow simulation in the present study takes cloud physics and precipitation processes into account.

Figure 15 shows the wind velocity vectors at the wind turbine height (60 m above the ground surface) in the smallest domain, d04 (Figure 10), at 9:40 a.m. JST on November 13, 2015. This figure confirms that the wind was easterly in the entire computational domain.

**Figure 14.** Computational domain used in the Weather Research and Forecasting (WRF) mesoscale model.

**Figure 15.** Distribution of the horizontal wind vectors in Domain 4, approximately 60 m above the ground surface. 9:40 a.m. JST, Nov. 13, 2015.

#### **5. Overview of the Numerical Simulation Method Based on Large-Eddy Simulation (RIAM-COMPACT)**

#### *5.1. Setting of Numerical Parameters*

As mentioned above, when easterly wind occurred, turbulence intensity was high around wind turbine #10, and thus the DEL of the wind blades was high from the viewpoint of actual measurement data analysis. Consequently, we inferred that the major cause was the influence of Mt. Benzaiten (elevation 519 m), located on the eastern side at a distance of approximately 300 m from wind turbine #10 (see Figure 16a). Thus, high-resolution numerical simulations for wind conditions were performed based on large-eddy simulation (LES) to scrutinize three-dimensional airflow structures, which directly affected the wind turbine blades due to terrain-induced turbulence. To compare flow distributions, the northerly wind was also examined. The northerly wind had the highest frequency of occurrence and was scarcely affected by terrain-induced turbulence (see Figure 16b).

For the numerical simulations, the RIAM-COMPACT (Research Institute for Applied Mechanics, Kyushu University, Computational Prediction of Airflow over Complex Terrain) natural terrain version was used, for which a collocated grid in a general curvilinear coordinate system was adopted. The collocated grid system is characterized by functions to define physical velocity components and pressure at the cell center of computational grids and to define variables at the cell faces, which are obtained with a contravariant velocity component multiplied by the determinant of the Jacobian. The numerical calculation method was based on the finite-difference method (FDM), and large-eddy simulation (LES) was employed to examine the turbulence models. The computational algorithm confirmed the fractional step method (F-S method) [25], and the time marching method was based on an explicit scheme of the Euler method. Poisson's equations regarding pressure were solved using successive over-relaxation (SOR). For the discretization of all the spatial terms, with the exception of the convective term in the equation, a second-order accurate central-difference scheme was applied. For the convective term, a third-order upwind difference scheme was applied. The interpolation technique was used for a fourth-order finite-difference that comprised the convective terms [26]. For the weighting of the numerical diffusion term in the convective term discretized by third-order upwind differencing, α = 0.5 is used as opposed to α = 3.0 [27] from the Kawamura-Kuwahara scheme in order to minimize the influence of numerical diffusion. The LES sub-grid-scale model employed the standard Smagorinsky model [28] combining a wall-damping function, setting the model coefficient as 0.1. In the present study, the LES is assumed to reproduce the wind tunnel testing. Therefore, the effects of atmospheric stability associated with vertical thermal stratification of the atmosphere and inflow turbulence were neglected. In addition, as in [12,16,17], the effects of surface roughness were taken into consideration by reconstructing surface irregularities in high resolution. A comparison between the Reynolds-averaged modeling (RANS) results and the present LES results is summarized in a recent article [13], and the prediction accuracy of the present LES approach by comparison with wind tunnel experiments is discussed in [19].

**Figure 16.** Comparison of the topographic section. (**a**) Northerly wind. (**b**) Easterly wind.

Conditions for numerical wind simulations were as follows: the computational domain held a space of 12.0 (x), 5.0 (y) and 2.5 (z) km for the streamwise direction, spanwise direction, and vertical direction, respectively (see Figure 17). A buffer zone was established in the upstream end of the computational domain, in which the terrain irregularities were reduced by 95% to form flat terrain. Similarly, a buffer zone was added to the downstream end of the domain. In the computational region, the maximum and minimum altitudes were 523.5 m and 0 m, respectively. This simulation utilized terrain elevation data with 10 m surface imagery provided by the Geospatial Information Authority of Japan (GSI). Generated grids, including added grids of the marginal area upstream and downstream of the computational region, were approximately eight million points in number, 496 (x) × 201 (y) × 81 (z) in three-dimensional coordinates. Horizontal grid resolutions in the vicinity of the wind turbine were 10 m in the x- and y-directions. The minimum resolution of the vertical grid was approximately 1.5 m above ground level to enable smooth drawing.

**(b)** 

**Figure 17.** Computational grids and domain. (**a**) Enlarged view. (**b**) Overall view.

The power law distribution was provided with N = 7.0, which is shown in Figure 18, regarding the inlet boundary conditions of numerical wind simulations. The present study focused on the impacts of terrain-induced turbulence for discussion, and thus did not examine fluctuations in the inlet airflow. A free-slip boundary condition was applied to the side and upper walls. A convective outflow condition was applied to the outlet boundary. A no-slip condition was imposed on the ground surface. Characteristics of reference scales in the present study are shown in Figure 19. In the figure, h represents di fferences in altitude in the computational region, Uin represents streamwise wind velocity of the inlet boundary at the maximum altitude point, and ν represents the coe fficient of kinematic viscosity. On the basis of the two types of reference scales, the dimensionless parameter, Re, is the Reynolds number (= Uin h/ν). For this simulation, Re = 104. The time step was specified as Δt = 2 × 10−<sup>3</sup> h/Uin. Furthermore, identical simulation conditions were applied to both northerly and easterly winds. In the present study, we used a vector-parallel supercomputer system named NEC SX-ACE. The NEC SX-ACE provides four cores per node. When using one node of this system, about 8 million grid points of LES simulation took several hours.

**Figure 18.** Inflow condition.

**Figure 19.** Two characteristic scales (Uin and h).

#### *5.2. Flow Visualization of Simulation Results*

Figure 20 illustrates the distribution of the streamwise (x) wind velocity component (u) in an instantaneous flow field. With the attentive observation of the two shadings in the figures, it is evident that the airflow pattern that was generated around the wind turbine in Figure 20b was significantly distinct in easterly winds compared to that of the northerly winds in Figure 20a Conclusively, with the visual e ffects of these illustrations, it was evident that when the wind blew from the east, wind turbine #10 was directly a ffected by the separated flow from the east, which was the terrain-induced turbulence formed due to Mt. Benzaiten (elevation 519 m), which is located upstream of the wind turbine.

**Figure 20.** Distribution of the streamwise wind velocity component on a vertical cross-section that includes wind turbine #10 and the instantaneous flow field. (**a**) Northerly wind. (**b**) Easterly wind.

#### *5.3. Proposal of Turbulence Evaluation Index (Uchida-Kawashima Scale\_1)*

The present study proposes two types of new indexes regarding wind conditions and loads. The former, the Uchida-Kawashima Scale\_1 (the U-K scale\_1), which is the turbulence evaluation index regarding wind conditions, is defined as formula (4) below; the latter, the Uchida-Kawashima Scale\_2 (the U-K scale\_2), which evaluates the fatigue damage, is defined as the fatigue damage evaluation index. Its definition is described later. The U-K scale\_1, which is the turbulence evaluation index, is obtained as follows: standard deviation is calculated using the streamwise (x) wind velocity component (u) at the hub height, and this is normalized by Uin, wind velocity at the maximum height above the ground of the inflow boundary, as shown in Figure 15. Instead of the average wind velocity at the hub height, the Uin inflow velocity at the inflow boundary is used for normalization; the formula is generalized and does not depend on wind directions or terrain undulations.

$$\text{U-K.Scale\\_1} = \frac{\sigma\_u}{\overline{\mathbf{U}\_{\text{in}}}} = \frac{\sqrt{\frac{1}{n}\sum\_{i=1}^{n}\left(\mu\_i - \overline{\mu}\right)^2}}{\mathbf{U}\_{\text{in}}} \tag{4}$$

#### *5.4. Analysis of Turbulence Statistics and Uchida-Kawashima Scale\_1 Verification Part 1*

Figure 21 shows time-series data (dimensionless time, 100) of the streamwise wind velocity component at the hub height (60 m above the ground) of wind turbine #10. The research data results of wind velocity revealed that the mean speed was lower, and the fluctuation was greater for easterly winds (shown in blue) than for northerly winds (shown in red). The threshold of the U-K scale\_1 was determined as follows: standard deviation, which was calculated with the streamwise (x) wind velocity component (u) assessed at the hub height, was normalized by inflow velocity at the maximum height above the ground of inflow boundary; the U-K scale\_1 is the turbulence evaluation index defined in formula (4). The value obtained for the northerly wind in this manner was 0.17, and similarly, that for the easterly wind was 0.25. Based on these data results, the present study defines the threshold of the U-K scale\_1 as 0.2. The validity of the value 0.2 of the threshold is discussed later.

Vertical profiles of turbulence statistics at the site of wind turbine #10 are shown regarding both northerly and easterly winds in Figures 22 and 23. Firstly, the mean values of the streamwise wind velocity components in Figure 22 were our focus. The profile of inflow wind velocity is also illustrated in green in Figure 22. For the northerly wind (in red), a speed-up effect of about 1.3 times (= 0.97/0.74) was obtained at the hub height (approximately 60 m) due to the topography effect on the wind. On the other hand, in the case of the easterly wind shown in blue, it was decelerated by about 0.4 times (= 0.29/0.74) under the influence of Mt. Benzaiten (elevation 519 m). Furthermore, it was confirmed in the figure that extremely large-scale velocity shears existed in the range of the swept area of the wind turbine (z\* = 20–100 m)(The variable z\* is the height above the ground). This large-scale velocity shear induces malfunctions of major parts that make up the wind turbine, such as the main shaft and gearbox. Thus, utmost caution must be exercised.

Secondly, among the three-dimensional components of standard deviations, the vertical profile is observed, as shown in Figure 23. The threshold value of 0.2 in the U-K scale\_1, which was previously described, is drawn in Figure 23. Results of the easterly wind, shown in blue, were greater than those of the northerly wind, shown in red, in all three components. Focusing on the wind turbine hub height (60 m above the ground), it is worth mentioning in particular that both northerly and easterly winds had almost the same level of values for the three components.

**Figure 21.** Time-series data of streamwise wind velocity from the numerical simulations. Red: northerly wind, Blue: easterly wind.

**Figure 22.** Vertical profiles of the streamwise wind velocity at wind turbine #10, with a time-averaged flow field. Red: northerly wind, Blue: easterly wind.

**Figure 23.** Vertical profiles of the non-dimensional standard deviations at wind turbine #10, with a time-averaged flow field. Red: northerly wind, Blue: easterly wind. (**a**) Streamwise direction. (**b**) Spanwise direction. (**c**) Vertical direction.

Compared with the threshold value of 0.2 in the U-K scale\_1, northerly patterns were lower than 0.2 in all three components; however, easterly patterns exceeded 0.2 in all three components. Based on these findings, it was quantitatively shown that the terrain-induced turbulence of the easterly wind had a three-dimensional structure, and turbulence was generated due to the presence of Mt. Benzaiten (elevation 519 m).

To verify the validity of the threshold value of 0.2 in the U-K scale\_1, configurations of inflow wind velocity profiles were altered using calculations of both northerly and easterly winds. Specifically, calculations using N = 4.0 and N = 10.0 were performed, and the obtained values were compared with the values of N = 7.0. These results are shown in Table 6. This table also shows the case of the threshold value of 0.2. The following results were obtained: in the case of the easterly wind, all of the calculated values were higher than 0.2 of the threshold values, and in all cases of the northerly wind, the calculated results were lower than 0.2 of the threshold value. It was suggested that the U-K scale\_1 did not depend on inflow wind velocity profiles, and thus 0.2 of the threshold level was a pertinent judgment criterion on the whole.

**Table 6.** Comparison of the values of the U-K Scale\_1 at wind turbine hub height (z\* 60 m) under different N values.

=


#### *5.5. Mt. Benzaiten Impact Assessment and Uchida-Kawashima Scale\_1 Verification Part 2*

The present study performed a computational simulation where the effects of Mt. Benzaiten (elevation 519 m) were eliminated from the examined impacts. Specifically, all the terrains higher than the elevation of wind turbine #10 were eliminated from the simulation. Furthermore, the validity of the threshold level in the U-K scale\_1 was considered. The grid resolution was changed from 10 m, which was used in previous numerical wind condition simulations in the horizontal direction, to five meters, which was half of the previous resolution. We investigated how this alternation in grid resolution affected the calculation results; that is, how they affected the values of the U-K scale\_1.

Figure 24 shows two types of situations for the easterly wind: simulation results of the current situation with Mt. Benzaiten (Figure 24a) and simulation results with Mt. Benzaiten eliminated (Figure 24b). In both simulations, the figures show the streamwise (x) wind velocity component (u) profiles and wind velocity vector in the vertical direction at the wind turbine #10 site for inlet wind to wind turbine #10 in the instantaneous flow field. From Figure 24a, for the current situation results, the separated flow of terrain-induced turbulence, which was generated from Mt. Benzaiten located in the upper stream of wind turbine #10, was precisely observed. Conclusively, wind velocity vectors at the wind turbine #10 site presented intricate distributions in height. Contrarily, the simulation results with Mt. Benzaiten removed (Figure 24b), as expected, showed that terrain-induced turbulence was not generated, and wind velocity vectors at the wind turbine site had an ideal distribution with a gradual increase in wind velocity by height.

(**a**) Simulation result of the current situation.

(**b**) Simulation result of removing Mt. Benzaiten (elevation: 519 m). 

**Figure 24.** Distribution of the streamwise wind velocity component on a vertical cross-section, which includes wind turbine #10 and wind velocity vectors at wind turbine #10, easterly wind, and instantaneous flow field.

Table 7 provides simulation results of the U-K scale\_1, where the horizontal grid resolution was changed to 5 m in high-resolution simulations to examine numerical wind conditions. The simulation with Mt. Benzaiten eliminated resulted in the U-K scale\_1 being 0.01, which was significantly lower than the threshold value of 0.2. In contrast, the simulation with the current situation, where terrain-induced turbulence exists, resulted in the U-K scale\_1 being 0.28, which was higher than 0.2 of the threshold

level, similar to the cases in Table 6. This meant that the effectiveness of the U-K scale\_1 and the validity of 0.2 of the threshold level were verified.


**Table 7.** Values of the U-K Scale\_1 with horizontal grid resolution set to 5 m. 

#### **6. In-situ Data Analysis of Impacts of Terrain-induced Turbulence on Fatigue Damage in Wind Turbine Blades**

#### *6.1. Relationship between Wind Speed and its Standard Deviation and Damage Equivalent Load (DEL)*

Figure 25 shows the results of the plotted nacelle wind speed and DEL obtained at 4 m/s or higher in wind turbine operation during the data measurement period of wind turbine #10 (3 November 2015, 0:00 a.m. JST–17 March 2016, 7:00 a.m. JST). In this figure, the values of the nacelle wind speed and DEL were in intervals of 10-min. The red symbol represents northerly wind, the blue symbol the easterly wind, and the black symbol represents design values calculated by the aero-elastic wind turbine simulation software Bladed. From these results, the following was clarified: firstly, when the speed of the easterly wind was between 6–10 m/s, the measured value exceeded the design value. This means that the blades of wind turbine #10 experienced wind loads that exceeded the design value due to terrain-induced turbulence, as a result of the proximity of Mt. Benzaiten (elevation 519 m), under 6–10 m/s easterly wind. Secondly, regarding the northerly wind, which had the highest frequency of occurrence during the above data measurement period, the damage was below the level of the design value in any wind speed class. Generally, it is known that wind speed and its standard deviation have a linear correlation. In Figure 26, the horizontal axis represents the standard deviation, which was converted by calculation from the wind speed axis of Figure 25. A strong correlation was revealed between DEL and the standard deviation values, as expected.

**Figure 25.** Relationship between wind speed (m/s) and damage equivalent load (DEL).

**Figure 26.** Relationship between standard deviation (m/s) and damage equivalent load (DEL). (**a**) Northerly wind. (**b**) Easterly wind.

#### *6.2. Investigation Regarding Accumulated Fatigue Damage of Wind Turbine Blades*

In the previous section, it was revealed that a linear trend was recognized between the nacelle wind speed (its standard deviation) and DEL in wind turbine operations at 4 m/s or higher, and it was possible to approximate to the regression line (see Figures 25 and 26). In other words, by assigning the DEL as a dependent variable (an objective variable) and nacelle wind speed as an independent variable (an explanatory variable), it is possible to apply a linear regression model between them. Figure 27 also shows a regression line of the nacelle wind speed and the DEL. The result of the northerly wind, which was represented by a red symbol, was defined as a low-turbulence flow case, and the results of the easterly wind, which was represented by a blue symbol, was defined as a high-turbulence flow case. At this point, the regression lines were modified at wind speeds of 10 m/s or higher. Similar to Figure 22, a black symbol was drawn as a design value calculated by the aero-elastic analysis software Bladed.

In the present study, as an index regarding the load, the Uchida-Kawashima Scale\_2 (the U-K scale\_2) was defined as the fatigue damage evaluation index, which was obtained by using two types of regression lines calculated based on the following actual measurement values and the design value on the basis of Bladed: one regression line was the northerly wind result as the low-turbulence flow case, and the other was the easterly wind result as the high-turbulence flow case. The definition of the U-K scale\_2 was as follows:

$$\text{U-K Scale\\_2} = \frac{\sum\_{i=1}^{n} DEL\_{\text{Proposal}}}{\sum\_{i=1}^{n} DEL\_{\text{Design}}} \tag{5}$$

At this point, we used m (S-N curve slope), which was used to calculate the DEL. The values raised to the *m*-th power were totaled, and the obtained value was raised to the (*1*/*m*)-th power for the integration of DEL, which is a scalar quantity. For example, the integration of two DELs is obtained with formula (6):

$$\text{DEL}\_{\text{total}} = \left(\text{DEL1}^m + \text{DEL2}^m\right)^{\frac{1}{m}} \tag{6}$$

Formula (5) refers to the ratio of the integrated value in the measured DEL to the integrated value in the designed DEL (Bladed). Therefore, the accumulated fatigue damage of wind turbine blades is as follows:

U-K scale\_2 > 1.0: more than the design value, impact of the terrain-induced turbulence: large. U-K scale\_2 - 1.0: design value and less, impact of terrain-induced turbulence: small.

**Figure 27.** Regression line between wind speed (m/s) and damage equivalent load (DEL).

In the present study, the data of both northerly and easterly winds were extracted from measured data for one year, April 2015–March 2016, corresponding to a 4 m/s or higher speed in the wind turbine's operation. Sequentially, using these extracted data and the U-K scale\_2, which was defined in formula (5), how terrain-induced turbulence influenced age-related degradation in wind turbine blades was assessed quantitatively. A total of 7485 data points (14.3%) were obtained for northerly wind and 2342 (4.5%) for easterly wind. In the case of northerly wind,

> U-K scale\_2 = 0.86 < 1.0: within the design value.

Contrarily, in the case of the easterly wind,

$$\text{U-K scale\\_2} = 1.60 > 1.0; \text{ over the design value.}$$

The U-K scale\_2 exceeded the design value, as was shown. Furthermore, the integrated value of fatigue damage was approximately 1.9 times greater for the easterly wind compared with the northerly wind. Based on this result, it was notably indicated that the blades of the wind turbine #10, which was the target turbine, were directly and intensely affected by the terrain-induced turbulence when easterly winds occurred.

#### **7. A Proposal for the Use of the U-K Scales and Future Research**

The present study showed that the effect of terrain-induced turbulence on wind turbine blades remains low as long as the value of the wind condition index, the U-K scale\_1, is equal to or smaller than the threshold value, 0.2. Furthermore, it was suggested that by combining this index and the threshold value of 1.0 or smaller for the U-K scale\_2, which is an index for wind loads on wind turbine blades, an optimal wind turbine siting can be planned with higher accuracy than with the conventional approach. Figure 28 shows an example flowchart for planning wind turbine siting with the use of the two generalized parameters proposed in the present study (the U-K scale\_1 and the U-K scale\_2). When planning the optimal wind turbine siting, it is desirable to maximize the output in prevailing winds, while simultaneously minimizing wind turbine failures. It is very likely that the use of the flowchart in Figure 28, even for a few prevailing wind directions, at a site under consideration will be effective. It is also possible to apply the flowchart for the placement of wind observation poles and for the so-called repowering of wind turbines at existing wind farms and individual wind turbine sites. Repowering refers to rebuilding old wind turbines in order to increase the wind power generation capacity and improve the efficiency of wind power generation.

The present study showed regression lines for the relationship between the DEL and wind speed at the nacelle height only for northerly winds as a low-turbulence flow case and easterly winds as a high-turbulence flow case. However, similar regression lines have also been evaluated for other wind directions. Thus, with the use of these regression lines, we are planning to quantitatively evaluate the effect of long-term DEL accumulation, for all wind directions, on the blades of wind turbine #10, which was the turbine studied in the present study, in the future. In addition, since WT #10 is affected by the wakes of other wind turbines in some wind directions, we are also planning to conduct analyses on the effect of wind turbine wakes on the blades of wind turbine #10.

**Figure 28.** An example of wind energy resource assessment based on the two reference scales (the U-K scales).

Our future studies will expand the present empirical research and aim to develop an advanced micro-siting technique. Specifically, this micro-siting technique will be able to 1) maximize the output from a group of wind turbines both onshore and offshore; and also 2) enable accurate assessments and predictions of fatigue damage and fatigue life of the main shaft, gearbox, and other main components of a wind turbine that result from terrain-induced turbulence and wind turbine wakes.
