Numerical Reproducibility of Terrain-induced Turbulence in Complex Terrain by Large-eddy Simulation ( LES ) Technique

In the present study, field observation wind data from the time of the wind turbine blade damage accident on Shiratakiyama Wind Farm were analyzed in detail. In parallel, high-resolution LES turbulence simulations were performed in order to examine the model’s ability to numerically reproduce terrain-induced turbulence. The comparison of the observed and simulated time series (1 second average values) from a 10 minute period from the time of the accident led to the conclusion that the settings of the horizontal grid resolution and time increment are important to numerically reproduce the terrain-induced turbulence that caused the wind turbine blade damage accident on Shiratakiyama Wind Farm. A spectral analysis of the same set of observed and simulated data revealed that the simulated data reproduced the energy cascade of the actual terrain-induced turbulence well.


Introduction
Recently, the number of accidents which involve wind turbines constructed on complex terrain in the mountains is increasing rapidly.Recent studies by the author of the present study and others 1), 2) indicate that the accidents they investigated are strongly associated with terrain-induced turbulence.In these studies, terrain-induced turbulence is defined as "temporal and spatial fluctuations of airflow that are mechanically generated due to terrain irregularities."Furthermore, the author of the present study classified terrain-induced turbulence roughly into two kinds.
The first kind is "extraordinary" terrain-induced turbulence which is generated with the passage of a typhoon or rarely occurring meteorological phenomena.That is, this kind of turbulence is terrain-induced turbulence which is commonly generated in wind directions that are different from the prevailing wind direction and occur infrequently throughout a year.It has been reported that this kind of terrain-induced turbulence caused a serious accident which involved cracks on wind turbine blades and other damage 1) .
The other kind of terrain-induced turbulence is "ordinary" terrain-induced turbulence which is generated under the prevailing wind direction.This turbulence has caused reduced power output from wind turbines and damage to the interior and exterior of wind turbines (e.g., breakdown of yaw motors and yaw gears), which have become evident issues 2) .
Uchida 1) investigated "extraordinary" terrain-induced turbulence.When Typhoon No. 0918 (Melor) passed the southern part of the main island of Japan between October 7 and 8, 2009, strong winds with extremely strong turbulence fluctuations were observed over the ridge of Mt.Shirataki and surrounding ridges in Houhoku-cho, Shimonoseki City, Yamaguchi Prefecture, Japan.This strong wind caused damage to a turbine blade on Shiratakiyama Wind Farm owned by Kinden Corporation (Fig. 1).Uchida 1) studied the cause of this accident by reconstructing the atmospheric phenomena occurring on spatial scales between a few meters and a few hundred kilometers using a computer simulation.In this study, the airflow field from the time of the accident was initially reconstructed in detail using a combination of a mesoscale meteorological model and RIAM-COMPACT, which is based on a large-eddy simulation (LES) turbulence model.Subsequently, the airflow fields in the vicinity of wind turbine blades were reconstructed separately using a RANS (Reynolds Averaged Navier-Stokes equation) turbulence model in order to evaluate the wind pressure on the wind turbine blades.For this analysis, the time-averaged flow field data from the LES simulation were used for the boundary conditions.Finally, stress exerted on the blades was calculated using the finite element method (FEM) with the RANS analysis results as the boundary conditions.

Fig.1 Photos of the blade damage
The results of the above analyses revealed that large values of stress occurred at the junction between the dorsal and ventral sides of the leading edge (LE) (the LE dorsal-ventral junction, hereafter) and that the locations at which these values occurred matched the locations of the actual damage.The static strength of the adhesive bond used for the LE dorsal-ventral junction is quite high.However, adhesive bonds will likely break when repeated loading as the one discussed above is applied to LE dorsal-ventral junctions as a result of low-cycle fatigue, that is fatigue associated with approximately 10 4 or fewer cycles for the number of cycles-to-failure.
During the passage of Typhoon 0918, a large number of wind direction deviations were observed under high wind conditions.Therefore, it was surmised that damage such as cracks formed at the LE dorsal-ventral junction as a result of low-cycle fatigue and this damage propagated further during subsequent wind turbine operation.
In contrast, Uchida 2) probed into "ordinary" terrain-induced turbulence by studying Taikoyama Wind Farm (located on Mt.Taikoyama, Nomura, Ine Town, Yosa District, Kyoto Prefecture, Japan).At this wind farm, at around 19:30 on March 12, 2013, a major accident occurred in which the nacelle and blades (hub height: approx.50.0 m; weight: approx.45.0 t) of Wind Turbine No. 3 fell to the ground as a result of a rupture of the tower in the vicinity of the tower top flange.Since the wind speed at the time of the accident was about 15.0 m/s and was within the design limit, investigations were conducted based on a standpoint that metal fatigue was the main cause of the accident.As a result of detailed investigations of the accident, the following was revealed: tensile stress on the welded joint between the tower shell and the tower top flange increased significantly because of the breakage of tower top bolts, which led to the formation of fatigue cracks on the internal wall of the tower shell in the vicinity of the weld toe.The formation of these fatigue cracks, in turn, caused the nacelle and blades to fall to the ground.Detailed numerical wind simulation results for the prevailing wind direction showed that many wind velocity shears which deviated from those predicted by a power law occurred at all of the wind turbine sites including Wind Turbine No. 3, the nacelle and blades of which fell to the ground.The simulation results also revealed that large wind direction deviations in the yaw direction of the wind turbines occurred frequently.Furthermore, the streamwise component of the wind velocity was relatively large, and the standard deviation of the vertical component of the wind velocity was large.The values of the standard deviation of the spanwise component of the wind velocity were approximately the same as those of the streamwise component of the wind velocity.From these findings, the following was conjectured.at a wind power generation facility, there is no doubt that the prediction and evaluation of turbulence intensity over complex terrain by CFD (computational fluid dynamics) such as LES and RANS will become increasingly important [3][4][5][6][7][8][9][10][11][12][13][14] . Gnerally, in RANS models, such as k -ε models, the standard deviation of the streamwise wind velocity which is attributable to terrain and/or surface roughness, Surf u  , is calculated using the values of turbulence kinetic energy, k; then, the standard deviation of the streamwise wind velocity, u  , is calculated by taking into account the background atmospheric turbulence intensity, where the background atmospheric turbulence intensity, a I , is taken as 0.1.The standard deviations of the spanwise and vertical wind velocities are calculated automatically using the ratios of these standard deviations to u  based on, for instance, an International Electrotechnical Commission (IEC) standard.However, this approach does not make it possible to evaluate the standard deviations of the three wind velocity components which result from turbulence structures that develop over complex terrain and thus deviate from those over flat, homogeneous terrain.In contrast, LES, which allows unsteady simulations, does not involve the above-mentioned assumptions and makes it possible to directly evaluate the standard deviations of the three wind velocity components using time series data of the three wind velocity components as is the case for evaluations of the standard deviations with field wind observation data.Therefore, LES is a highly effective means to predict and evaluate turbulence intensity over complex terrain.
In the present study, field observation wind data from the time of the wind turbine blade damage accident at Shiratakiyama Wind Farm, which was investigated in Uchida 1) , are analyzed in detail.In addition, a high-resolution LES turbulence simulation is performed by refocusing attention on the connection between the blade damage and the terrain-induced turbulence which caused the blade damage.Based on the results from this simulation, the numerical reproducibility of terrain-induced turbulence by the LES model is examined.

Shiratakiyama Wind Farm and Airflow Characteristics from the Time of the Blade Damage Accident
The Shiratakiyama Wind Farm owned by Kinden Corporation is located on a ridge near Mt.Shirataki in Houhoku-cho, Shimonoseki City, Yamaguchi Prefecture, Japan.On this wind farm, twenty 2,500 kW wind turbines manufactured by General Electric Company (wind turbine hub height: 85.0 m; swept-area diameter: 88.0 m; height of the upper end of the swept area above the ground surface: 129.0 m) are deployed (see Figs. Peer-reviewed version available at Energies 2018, 11, 2638; doi:10.3390/en11102638 Fig. 5 shows the trajectory of Typhoon No. 0918 (Melor), pressure charts from 9:00, October 7, 2009 and 9:00, October 8, 2009, and the number of wind turbine incidents on these two days.Generally, strong winds blow to the east of the direction of motion of a typhoon.However, on October 7, 2009, on which the blade damage accident occurred, the center of the typhoon was located to the south of the main island of Japan as shown by the trajectory of the typhoon.Therefore, the position of Shiratakiyama Wind Farm on this day was in the range between north and west of the center of the typhoon.At that time, a steep horizontal pressure gradient (i.e., a zone with densely packed isobars extending from east to west) was present between the typhoon and the high pressure located to the north of the typhoon (Fig. 5).As a result, the wind was flowing from the east around the typhoon to the north of the typhoon, which caused strong, north-easterly wind to flow into the site of the blade damage accident, causing the accident.In the present study, a new high-resolution LES turbulence simulation is performed to investigate numerical reproducibility of the terrain-induced turbulence which caused the wind turbine blade damage accident.

Overview of Numerical Simulation Method
The present study uses the RIAM-COMPACT natural terrain version software package 1,2,[15][16][17][18][19][20][21][22][23][24] , which is based on a collocated grid in a general curvilinear coordinate system. I this collocated grid, the velocity components and pressure are defined at the grid cell centers, and variables that result from multiplying the contravariant velocity components by the Jacobian are defined at the cell faces.For the numerical technique, the finite difference method (FDM) is adopted, and a large-eddy simulation (LES) model is used for the turbulence model.
In the LES model, a spatial filter is applied to the flow field to separate eddies of various scales into grid-scale (GS) components, which are larger than the computational grid cells, and sub-grid scale (SGS) components, which are smaller than the computational grid cells.Large-scale eddies, i.e., the GS components of turbulence eddies, are directly numerically simulated without the use of a physically simplified model.In contrast, dissipation of energy, which is the main effect of small-scale eddies, i.e., the SGS components, is modeled according to a physics-based analysis of the SGS stress.
For the governing equations of the flow, a filtered continuity equation for incompressible fluid (Eq.( 3)) and a filtered Navier-Stokes equation (Eq.( 4)) are used.In the present study, the effects of atmospheric stability associated with vertical thermal stratification of the atmosphere are neglected.As in Uchida 1,2) , the effects of the surface roughness are taken into consideration by reconstructing surface irregularities in high resolution.For the computational algorithm, a method similar to a fractional step (FS) method 25) is used, and a time marching method based on the Euler explicit method is adopted.The Poisson's equation for pressure is solved by the successive over-relaxation (SOR) method.For discretization of all the spatial terms except for the convective term in Eq. ( 4), a second-order central difference scheme is applied.For the convective term, a third-order upwind difference scheme is applied.The interpolation technique by Kajishima et al. 26) is used for the fourth-order central differencing that appears in the discretized form of the convective term.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 from the Kawamura-Kuwahara scheme 27) in order to minimize the influence of numerical diffusion.For LES subgrid-scale modeling, the standard Smagorinsky model 28) is adopted with a model coefficient of 0.

Overview of Numerical Simulation Set-up
The computational domain used in the present study extends over the space of 21.0 km (x) × 10.0 km (y) × 3.4 km (z), where x, y, and z are the streamwise, spanwise, and vertical directions, respectively 9).The maximum terrain elevation within the computational domain is 686.0 m, and the minimum terrain elevation is the sea surface (0.0 m).Terrain elevation data with 10.0 m and 50.0 m spatial resolutions from the Geospatial Information Authority of Japan are used.The total number of computational grid points, 421(x) × 201(y) × 51(z), is approximately 4.3 million, and the grid spacing is uniform (50.0 m) in both the x-and y-directions.Fig. 10 shows a comparison of terrain configurations constructed from two data sets, i.e., 10.0 m and 50.0 m resolutions, for the area of Shiratakiyama Wind Farm.The comparison reveals no significant difference, and the accuracy of the reconstructed terrain is approximately the same between the two data sets.The grid spacing is non-uniform in the z-direction so that the density of grid points increases smoothly toward the ground surface.The minimum vertical grid spacing is 2.0 m.An earlier study of the wind blade damage accident 1) conjectured that the accident occurred with northeasterly wind.Therefore, the simulations in the present study are also performed for northeasterly wind conditions.Regarding the boundary conditions, the wind velocity profile for surface roughness classification III is given as the inflow wind velocity at the inflow boundary (Fig. 11).At the side and upper boundaries, slip conditions are applied, and convective outflow conditions are applied at the outflow boundary.On the ground surface, a no-slip boundary condition is imposed.The non-dimensional parameter Re in Eq. ( 4) is the Reynolds number (= Uh/ν).For the simulations, Re = 10 4 is used.Fig. 12 illustrates the characteristic scales used in the present simulations: h is the difference between the minimum and maximum terrain elevations within the computational domain (= 686.0 m), U is the wind velocity at the inflow boundary at the height of the maximum terrain elevation within the computational domain, and ν is the kinematic viscosity of air.The time increment is set to Δt = 2 × 10

Simulation Results and Discussions
First, a comparison is made on the results from the simulations which used terrain elevation data with 10.0 m and 50.0 m resolutions from the Geospatial Information Authority of Japan (Section 4).Fig. 13 shows the vertical profiles of the mean streamwise wind velocity from the two simulations.Specifically, the plotted values were acquired by time-averaging (frame-averaging) the streamwise wind velocity at Wind Turbine No. 17 over the non-dimensional time period of t = 200.0-400.0.These values correspond to output from a RANS model.The variable z* on the vertical axis represents the height above the terrain surface (m), and the horizontal axis represents the averaged streamwise wind velocity normalized by the inflow wind velocity U (m/s).As can be presumed from the fact that no significant differences exist between the two terrain data sets (Fig. 10), the tendencies of the results from the two simulations including those plotted in Fig. 13 were nearly the same.In light of this finding, further discussions are presented based on the results from the simulation which used the 10.0 m resolution terrain elevation data.flow field, in which turbulent eddies are distorted significantly away from isotopic forms.This finding is also consistent with the large number of wind direction deviations which were recorded in the field observation wind data.Therefore, as already mentioned in the introduction, it is quite difficult to accurately evaluate values of the standard deviations the three wind velocity components with the use of RANS models such as k -ε models for turbulent airflow that develops over complex terrain.Specifically, with this approach, the standard deviation of the x-component of the wind velocity is determined from turbulence kinetic energy k, and the standard deviations for the y-and z-components of the wind velocity are automatically calculated using the ratios of these standard deviations to the standard deviation of the x-component of the wind velocity (see Section 1).For the field observation wind data set, scalar horizontal wind speed, Uscalar, was measured and recorded.This scalar horizontal wind speed, Uscalar, can be related to the horizontal wind speed acquired from the simulation as Uscalar = (u 2 +v 2 ) 1/2 , where u and v are the streamwise and spanwise wind velocities acquired from the simulation.
Since the mean horizontal wind speed acquired from the field wind observations was 8.40 m/s and the mean horizontal wind speed from the numerical simulation was 0.27, all of the wind velocity data values from the numerical simulation were multiplied by 8.40/0.27 to ensure that the mean horizontal wind speed from the numerical simulation would be equal to 8.40.Accordingly, the inflow wind velocity was set to approximately 31.0 m/s for the numerical simulation.This value is in close agreement with the value of the wind velocity obtained from the mesoscale meteorological model in Uchida   Peer-reviewed version available at Energies 2018, 11, 2638; doi:10.3390/en11102638

Conclusion
In the present study, field observation wind data from the time of the wind turbine blade damage accident on Shiratakiyama Wind Farm, which was studied by Uchida 1) , were analyzed in detail.In parallel, high-resolution LES turbulence simulations were performed in order to examine the model's ability to numerically reproduce terrain-induced turbulence.First, a comparison was made between terrain elevation data with 10.0 m and 50.0 m spatial resolutions from the Geospatial Information Authority of Japan.When a uniform grid spacing of 50.0 m was used for the horizontal grid for the airflow calculation, the accuracy of the simulation results was approximately the same between the simulations which used the two data sets.As a result, no significant difference identified between the vertical profiles of mean streamwise wind velocity simulated over the terrains constructed from the two terrain data sets.Subsequently, vertical profiles of the standard deviations of the three wind velocity components at Wind Turbine No. 17 were examined.This examination revealed that the value of the standard deviation for the y-component of wind velocity exceeded that for the x-component of wind velocity at the wind turbine hub height (= 85.0 m) and below.This finding is consistent with the large number of wind direction deviations which were recorded in the field observation wind data.The comparison of the observed and simulated time series (1 second average values) from a 10 minute period from the time of the accident led to the conclusion that the settings of the horizontal grid resolution and time increment are important to numerically reproduce the terrain-induced turbulence that caused the wind turbine blade damage accident on Shiratakiyama Wind Farm.A spectral analysis of the same set of observed and simulated data revealed that the simulated data reproduced the energy cascade of the actual terrain-induced turbulence well.

Fig. 2
Fig.2 Location of Shiratakiyama Wind Farm in Yamaguchi Prefecture, Japan

Fig. 7 Fig. 7 Fig. 8
Fig. 7 shows time series data from one of the wind vanes and one of the anemometers mounted on the nacelle of Wind Turbine No. 17 (50.0m above the ground surface).In Figs 7(a) to 7(d), the horizontal axes indicate time (Midnight, October 7, 2009 to noon, October 8, 2009).These figures also include indications of 1) the time period in which the wind turbine blade damage accident occurred (14:00 -22:00, October 7, 2009) and 2) the start time of the tower vibration (15:50).On October 7, 2009, the wind speed started increasing gradually in the area of Shiratakiyama Wind Farm in the afternoon.As described earlier, at around 14:00, October 7, 2009, the time at which the wind turbine blade damage accident began, north-easterly wind was flowing into the farm near the

Fig. 9
Fig.9 Computational domain and grid

Fig. 14 shows
Fig.14shows the vertical profiles of the standard deviations of the three wind velocity components at Wind Turbine No. 17.The standard deviations were calculated with respect to the time-averaged (frame-averaged) wind velocity components from the non-dimensional time period of t = 200.0-400.0.In the present study, the fluctuating wind velocity components (gusts) that are present in the observed inflow wind are not included in the inflow wind in the simulation, thus only the airflow fluctuations of terrain-induced turbulence which are generated due to the terrain irregularities are evaluated.The left panel of Fig.14shows the entire range, and the right panel of Fig.14shows an enlarged view for the range of z* = 0.0 -200.0 m.The left panel (entire range) indicates that the maximum value of the standard deviation of the x-component of the wind velocity occurs in the vicinity of z* = 300.0m.The right panel (enlarged view) shows that the values of the standard deviations of the three wind velocity components are relatively large within the swept area of the wind turbine.In particular, it should be especially mentioned that the values for the y-component exceed those for the x-component at the wind turbine hub height (= 85.0 m) and below.This finding suggests the formation of an anisotropic turbulent

Fig 13 .
Fig 13.Profiles of mean streamwise wind velocity at Wind Turbine No. 17

Fig. 14
Fig.14 Profiles of the standard deviations of the three wind velocity components at Wind Turbine No.17, resolution of terrain elevation data = 10.0 m 1) .The non-dimensional time on the horizontal axis t (= 200.0 -300.0) was converted to full scale (s) based on T = t (h/U), where h = 686.0m and U = 31.0m/s for this case.As a result, the non-dimensional time interval of 100 was converted to approximately 2,209 seconds (approximately 37 minutes) in full scale (time step interval: 0.044 s).In Fig.15, the red solid line indicates field wind observation data (scalar horizontal wind speed) for 600 seconds (10 minutes), and the black solid line indicates the horizontal wind speed calculated from the numerical simulation for approximately 2,209 seconds (approximately 37 minutes).An examination of the numerical simulation results reveals the following.Rapidly fluctuating high-frequency components which are present in the field wind observation data are not fully reproduced in the simulation data.However, wind velocity fluctuations with very short periodic cycles which are attributable to generation of terrain-induced turbulence are reproduced for the most part.As a consequence, both the standard deviation of the horizontal wind speed (m/s) and turbulence intensity evaluated from the field observation and simulated wind data are successfully in close agreement (Fig.15).In order to obtain the time series with the above-mentioned periodicity, settings of parameters Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 4 July 2018 doi:10.20944/preprints201807.0077.v1Peer-reviewed version available at Energies 2018, 11, 2638; doi:10.3390/en11102638that are discussed below (i.e., horizontal grid resolution and time increment) are particularly important.When strong north-easterly wind flows into Shiratakiyama Wind Farm, airflow characterized by rapid fluctuations in wind direction and speed (terrain-induced turbulence), which originates from a relatively large terrain feature located upwind of the wind farm, flows into the wind turbines (see Uchida1) ).Therefore, it is necessary to accurately reproduce 1) a relatively large-scale terrain feature which is present upwind of the wind farm and 2) the three-dimensional structure of terrain-induced turbulence which is generated due to the terrain feature by sufficiently resolving the turbulence both temporarily and spatially.To achieve this goal, a grid resolution of approximately 10.0 m is required for horizontal cross sections.in order to capture unsteady fluid properties of terrain-induced turbulence, a sufficiently small time increment (Δt = 2×10 -3 h/U in the present study) is required.To investigate how accurately the terrain-induced turbulence simulated in the present study has reproduced the terrain-induced turbulence over the wind farm in frequency space, a power spectral analysis is performed on the fluctuating components of the time series in Fig. 15 using a fast Fourier transform (FFT).The results of this analysis are shown in Fig. 16.In this figure, the vertical axis shows the power spectra, which is non-dimensionalized by frequency f (Hz) and standard deviation σ (m/s), and the horizontal axis shows frequency f (Hz).The power spectra show that the simulated data have reconstructed the energy cascade of the actual terrain-induced turbulence over the principal wavenumber band well.

Fig. 15
Fig. 15 Comparison of horizontal wind speed between observed data (red line) and simulated data converted to full scale (black line)

Finally, using the
dominant frequency from Fig.16, f = 0.04 Hz, the elevation of Tenjogadake, h = 691.1 m (Fig.10), and the choice of the inflow wind velocity U = 31.0m/s together yields 0.89 for the non-dimensional frequency, St (= f h/U).This value nearly agrees with the value of the vortex shedding frequency, St = 0.87, which was obtained from wind tunnel experiments with simple topography (a two-dimensional ridge and a three-dimensional isolated hill)24) .Thus, it is likely that the terrain-induced turbulence which caused the wind turbine blade damage accident on Shiratakiyama Wind Farm was attributable to rapid wind speed and direction fluctuations which were caused by vortex shedding from Tenjogadake (elevation: 691.1m) located upstream of the wind farm.Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 4 July 2018 doi:10.20944/preprints201807.0077.v1

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 4 July 2018 doi:10.20944/preprints201807.0077.v1 Peer
The exciting force on Wind Turbine No. 3 increased due to the -reviewed version available at Energies 2018, 11, 2638; doi:10.3390/en11102638effect of terrain-induced turbulence.As a result of the increased exciting force, additional load was imposed in the vicinity of the tower top flange, and thus increased metal fatigue in multiple bolts.In light of the recent increase in wind turbine accidents such as the ones described above, laws and regulations about wind power generation in Japan have been reviewed and amended.Specifically, the Nuclear and Industrial