High-Resolution LES of Terrain-Induced Turbulence around Wind Turbine Generators by Using Turbulent Inflow Boundary Conditions

We have developed an LES (Large-Eddy Simulation) code called RIAMCOMPACT (Research Institute for Applied Mechanics, Kyushu University, Computational Prediction of Airflow over Complex Terrain). The analysis domain of this numerical model extends from several meters to several kilometers. The model is able to predict airflow over complex terrain with high accuracy and is also now able to estimate the annual power output of wind turbine generators with the use of field observation data. In the present study, a numerical simulation of turbulent airflow over an existing wind farm was performed using RIAM-COMPACT and high-resolution elevation data. Based on the simulation results, suitable and unsuitable locations for the operation of WTGs (Wind Turbine Generators) were identified. The latter location was subject to the influence of turbulence induced by small topographical variations just upwind of the WTG location.


Introduction
We have developed an unsteady and non-linear wind synopsis simulator called RIAM-COMPACT (Research Institute for Applied Mechanics, Kyushu University, Computational Prediction of Airflow over Complex Terrain) in order to simulate the airflow on a microscale, i.e., a few tens of km or less [1]- [10].In RIAM-COMPACT, large-eddy simulation (LES) has been adopted for turbu-T.Uchida DOI: 10.4236/ojfd.2017.74035512 Open Journal of Fluid Dynamics lence modeling.LES is a technique in which the structures of relatively large eddies are directly simulated and smaller eddies are modeled using a sub-grid scale model.Other wind-synopsis software which has been developed and is currently available within Japan utilizes models based on Reyolds-averaged Navier-Stokes equations because of the computational time requirement.These models are called RANS models and are utilized for simulations of stationary flow fields in which the flow properties remain constant in time [11].However, computer performance has improved rapidly in recent years, and the constraints imposed by computational time have been drastically reduced.LES models, which solve spatially-averaged Navier-Stokes equations, are able to simulate unsteady wind fields subject to constant changes in flow properties.This aspect of LES models makes them considerably different from RANS models, in which the turbulent flow is temporally (Reynold's) averaged.Furthermore, the number of model parameters to be tuned in LES models is significantly smaller than that in RANS models, thus, LES models excel in versatility.If the characteristics of unsteady wind conditions can be easily predicted by an LES model and the simulation results can be captured visually by an animation or other means, the LES model can serve as an alternative to a costly wind tunnel experiment and can also contribute to the series of investigations conducted prior to the construction of a wind farm (WF).
Recently, the values of the availability factors of wind turbine generators (WTGs) on large WFs constructed on complex terrain have fallen below the values projected originally.In other words, WTGs with notably low power output have been identified, and the issues of both internal and external breakdowns of WTGs have emerged.The major cause of these problems is that slight topographical variations become a source of turbulence (terrain-induced turbulence), and turbulence is generated mechanically (directly) at these topographical variations.Furthermore, with global warming, interannual change in the average wind speed as well as the prevailing wind direction on the WF may become causes of the above-mentioned problems reported at large WFs.In our earlier studies, we performed detailed wind synopsis simulations for existing WFs with the use of elevation data with high resolution, i.e., less than 10 m.
Based on the simulation results, we proposed techniques which address issues concerning the relocation of WTGs within a WF such as if the current deployment locations of WTGs should be maintained and if poorly operating WTGs should be relocated or discontinued because of the operating cost.However, our past studies addressed only the influence of the topographical variations on the WTGs and their surroundings and curtailed investigations on the influences of turbulence which is generated upstream of the WF and flows into the WF and of surface roughness on WTGs and their surroundings.In the present study, in order to simulate the conditions in reality more accurately than in the previous studies, roughness blocks are placed upwind of a WF within the computational domain.This arrangement allows the WTGs under investigation to be com-    Therefore, the accuracy of these maps and data is significantly higher than that of, for instance, the 50 m elevation data of the Geographical Institute of Japan.Elevation data created by GIS such as those in Figure 2(a) can be readily used by RIAM-COMPACT.The time required for the construction of elevation data from a single paper map is only a few days.If information related to WTGs such as the actual dimensions of the rotor diameter and the hub-height, the wind directions at the swept areas, display colors, and the latitude and longitude of the locations of WTGs in decimal degrees based on the World Geodetic System are specified within RIAM-COMPACT, diagrams of WTGs can be inserted in the computational domain (see Figure 2).The rotor diameter and hub-height of the WTGs considered in the present study are 52 m and 44 m, respectively.For convenience, these WTGs will be referred as WTG-A and WTG-B and displayed with the letters A and B in figures hereafter.In order to achieve simulation results with accuracy higher than that from earlier simulations, a wind synopsis simulation is performed with RIAM-COMPACT with the use of high-resolution elevation data created with the above-mentioned method as input data.

LES Turbulence Simulation by RIAM-COMPACT
To simulate the airflow over complex terrain with high accuracy and avoid numerical instability, the simulation is performed with RIAM-COMPACT in which collocated grids in a general curvilinear coordinate have been adopted.In the collocated grid system, the velocity components and pressure are defined at the cell centers, and variables which result from the covariant velocity components multiplied by the Jacobian are defined at the cell faces.As for the computational technique, the finite-difference method is adopted, and LES is used for the turbulence model.In LES, a spatial filter is applied on the flow field to separate eddies of various scales into grid scale (GS) components, which are larger than the computational grids, and subgrid scale (SGS) components, which are smaller than the computational grids.Large-scale eddies, i.e. the GS components of turbulence eddies, are numerically simulated directly without relying on the use of a physically simplified model.The main effect of small-scale eddies, i.e.
the SGS components, is to dissipate energy, and this dissipation is modeled based on the physical considerations of the SGS stress.For the governing equations of the flow, a spatially-filtered continuity equation for incompressible fluid and a spatially filtered version of the Navier-Stokes Equation are used.Because the present study is on airflow prediction in high wind conditions, the effect of the temperature stratification which is generally present in the atmosphere is neglected.The computational algorithm and the time marching method are based on a fractional-step (FS) method [12] and the Euler explicit method, respectively.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, a second-order central difference scheme is applied.For the convective term, a third-order upwind difference scheme is applied.Generally, the third-order upwind difference equation consists of a fourth-order central difference term and a numerical diffusion term which takes the form of a fourth-order derivative.For the fourth-order central differencing, an interpolation technique based on 4-point differencing and 4-point interpolation by Kajishima [13] is used in the present study.In the weighting of the numerical diffusion term of the third-order upwind differencing, α = 3.0 is commonly applied in the Kawamura-Kuwahara Scheme [14].However, α is set to 0.5 in the present study to minimize the influence of numerical diffusion.For LES SGS modeling, the commonly used Smagorinsky model [15] is adopted.A wall-damping function is also used with a model coefficient of 0.1.
Regarding the boundary conditions, a slip condition is applied at both the upper and side boundaries.Specifically, for the upper boundary, the gradient of the horizontal wind speed components (u, v) in the vertical direction (z) and the vertical wind velocity component (w) are all set to zero.For the side boundaries, the gradients of the streamwise and vertical wind velocity components (u, w) in the spanwise direction (y) and the spanwise wind velocity component (v) are all set to zero.For the outflow boundary and the ground, a convective outflow condition and a non-slip condition are imposed, respectively.At the inflow boundary, a vertical profile of the horizontal streamwise wind speed, U in , is given using a 1/7 power law.

LES Turbulence Simulation by RIAM-COMPACT
Figure 3 shows the flow field over the entire computational domain, displaying the instantaneous wind velocity in the streamwise direction (x), i.e., u , at non-dimensional time t* = 100.Figure 4 shows the vertical profile of the mean streamwise wind velocity, U (= < u >), and those of the standard deviations of the individual wind velocity components, u  Figure 1.These figures suggest that separated flow occurs at the roughness blocks placed near the inflow boundary, advances downstream, and develops spatially.As a result, in the simulation, the WTGs are completely immersed within a turbulent boundary layer.Figure 5 shows vertical profiles of horizontal wind velocity vectors at and in the vicinity of the WTGs from non-dimensional time t* = 100.At this time, the flow is fully developed.Neither flow separation nor its accompanying vortex formation occurs upwind or downwind of WTG-A (Figure 5(a)).As for the vertical distribution of the horizontal wind velocity at heights across the swept area, the distribution is nearly uniform for WTG-A (Figure 5(b)).Figure 5(b) also indicates that the horizontal wind velocity is enhanced locally at some heights within the swept area of WTG-A due to topographical effects.In contrast, at WTG-B, the presence of a reverse flow and its accompanying vortex formation are evident beneath the swept area (indicated by the rectangle in Figure 5(a)).Large wind velocity deficits are also present at WTG-B as indicated by the arrow in Figure 5(b), and these deficits are attributable to the effects of the turbulence induced immediately upwind of WTG-B by the topographical variations (terrain-induced turbulence).More specifically, WTG-B is subject to the Subsequently, the simulated wind field is evaluated more quantitatively than in the previous section.Figure 6 shows the temporal change at WTG-A and WTG-B of the non-dimensional fluctuating component of the horizontal streamwise (x) wind velocity i.e., u' (= u U − ), where u and U are the instantaneous and the mean streamwise horizontal wind velocities at the 44 m hub-height, respectively.The velocity fluctuations at WTG-A in Figure 6(a) indicate the influence of turbulence in the wind flowing into the WF, suggesting that turbulence generated directly by the topographical variations (terrain-induced turbulence) is quite small.On the other hand, velocity fluctuations with significantly large amplitudes appear periodically at WTG-B (Figure 6(b)).These fluctuations correspond to the velocity fluctuations (turbulent fluctuations) associated with vortex shedding which occur periodically at the topographical undulation in the proximity of the WTG as discussed earlier.
Figure 7 and Figure 8 illustrate the vertical profiles of the mean streamwise wind velocity, U (= < u >), and the standard deviations of the streamwise, spanwise, and vertical wind velocity components, i.e, σ u ( 2  The power curve of a WTG is calculated using the value of the inflow velocity at the hub-center and without taking into account the presence of the WTG.
Furthermore, the calculation assumes that the values of the velocity shear present in the vertical wind profile are equal to those in a wind profile power law with exponents of approximately 5 to 7. Therefore, when the values of the vertical wind shear deviate significantly from those predicted by the wind profile power law, it can be anticipated that the actual output of the WTG falls significantly below the output predicted by the numerical model used to create the power curve.It is likely that very large velocity shear will be an important research topic in the future in connection with the issues of WTG tower vibration and the fatigue strength of WTGs.
The vertical profiles of the standard deviations of the three wind velocity components reveal that the values of the standard deviations at WTG-A (Figure 8(a)) are roughly the same as those for the wind flowing into the WF evaluated over flat terrain (Figure 4).Therefore, WTG-A is nearly free from the influence of terrain-induced turbulence.In contrast, at WTG-B, the maximum value of the   This series of fluid phenomena, which originates from the topographical variations, leads to mechanical (direct) generation of turbulence (terrain-induced turbulence) as shown in Figure 9.

Summary
In the present study, in cooperation with Eurus Energy Holdings Corporation, elevation data for a wind farm (WF) within Japan were created with high resolution, i.e., less than 10 m, and including the latest land development status.With these elevation data, a wind synopsis simulation was performed by taking into account the turbulence in the wind flowing into the WF.As found in our past studies, small topographical variations became a source of turbulence, and turbulence (terrain-induced turbulence) was mechanically (directly) generated at these topographical variations.In addition, suitable and unsuitable locations for WTGs were identified.Quantitative evaluation of the threshold for velocity fluctuations, above which adverse effects on the WTGs are anticipated, is a major research topic for the future.Advancement in this research topic will likely contribute significantly to 1) formulation of criteria for WTG siting risks, 2) construction of maps of WTG siting risks according to the formulated criterion, and 3) establishment of standards for WTG construction and regulations for WTG operation, which are to be used by the Japanese wind power generation industry (J-Class wind power generation).
It is anticipated that many future wind power generation facilities will need to be constructed in mountain areas or other locations which are less than ideal.
Accordingly, future evaluations of wind power generation facility projects require higher precision and accuracy than the evaluations of such projects conducted thus far.The use of the WF wind synopsis diagnosis (detailed airflow analysis) proposed in the present study which utilizes elevation data with high resolution, i.e., less than 10 m, is recommended for wind synopsis investigations conducted prior to WTG construction.In terms of the visualization of the wind field, wind velocity distributions on the hub-height horizontal cross-section have conventionally been evaluated.However, such evaluations alone will be insufficient for future evaluations of wind power generation facility projects and visualizations of the fluctuating wind velocities, standard deviations of the wind velocities, and vertical cross-sections of the horizontal wind speed and direction at

Figure 1
Figure1illustrates the computational domain used in the present study.The area in the computational domain represents an actual WF within Japan.The simulation is performed for the case of westerly wind.The dimensions of the computational domain are 6 km, 2 km, and 1.44 km in the streamwise (x), spanwise (y), and vertical (z) directions, respectively.The numbers of grid cells are 301, 101, and 51 in the x, y, and z directions, respectively.The horizontal grids are equally-spaced at Δx = Δy = 20 m.The vertical grids are unequally-spaced in the range of Δz = 0.85 -288 m.The characteristic length scale H, which is set to the difference between the minimum and maximum ground surface elevations within the computational domain, is 288 m in the present simulation.In order to simulate the turbulent boundary layer over the WF, five roughness blocks of height H (=288 m) are placed upstream of the WF in the computational domain.Figure 2(a) and Figure 2(b) show a plane view and a bird's-eye view, respectively, of the topography in the area with the WTGs.In the present study, the topography around the locations of the WTGs are reconstructed in detail with 5 m resolution from CAD data in the DXF format, and the 5 m resolution topographic data are joined smoothly to the 50 m elevation data of the surrounding area which are available from the Geographical Survey Institute of Japan (Figure 2(a)).The application of geographical information system (GIS) techniques

Figure 2 (Figure 1 .
Figure1illustrates the computational domain used in the present study.The area in the computational domain represents an actual WF within Japan.The simulation is performed for the case of westerly wind.The dimensions of the computational domain are 6 km, 2 km, and 1.44 km in the streamwise (x), spanwise (y), and vertical (z) directions, respectively.The numbers of grid cells are 301, 101, and 51 in the x, y, and z directions, respectively.The horizontal grids are equally-spaced at Δx = Δy = 20 m.The vertical grids are unequally-spaced in the range of Δz = 0.85 -288 m.The characteristic length scale H, which is set to the difference between the minimum and maximum ground surface elevations within the computational domain, is 288 m in the present simulation.In order to simulate the turbulent boundary layer over the WF, five roughness blocks of height H (=288 m) are placed upstream of the WF in the computational domain.Figure 2(a) and Figure 2(b) show a plane view and a bird's-eye view, respectively, of the topography in the area with the WTGs.In the present study, the topography around the locations of the WTGs are reconstructed in detail with 5 m resolution from CAD data in the DXF format, and the 5 m resolution topographic data are joined smoothly to the 50 m elevation data of the surrounding area which are available from the Geographical Survey Institute of Japan (Figure 2(a)).The application of geographical information system (GIS) techniques

and w σ ( 2 wσFigure 3 .
Figure 3. Contour plots of the streamwise (x) wind velocity component over the entire computational domain.The wind velocity range of −1.0 to 1.5 (Non-dimensional) is divided into 30 equal intervals.(Instantaneous field; Non-dimensional time t* = 100).(a) Side view, vertical cross-section through WTG-B.(b) View from aloft, horizontal cross-section at the hub-height (44 m above the ground surface).

Figure 4 .
Figure 4. Characteristics of the wind flowing into the WF evaluated at the location marked by the circle in Figure 1.The values plotted on the horizontal axes are normalized by the wind velocity, U in , at H = 288 m above the ground surface at the inflow boundary.The characteristic length scale, H = 288 m, adopted for the present study is illustrated in the figure.The vertical axes represent the height above the ground surface in full scale, z* (m).The values were computed over a period of non-dimensional time of 100.(a) Mean horizontal wind velocity profile.(b) Standard deviation of the individual wind velocity components.

Figure 5 .
Figure 5. Wind velocity vectors at and in the vicinity of the WTGs (Instantaneous field; Non-dimensional time t* = 100).(a) Vertical cross-sections of horizontal wind velocity at and in the vicinity of the WTG locations.(b) Vertical profiles of horizontal wind velocity at the WTG locations from the same time as in Figure 5(a).

Figure 6 .
Figure 6.Temporal change in the non-dimensional fluctuating horizontal streamwise (x) wind velocity, u', at the WTG hub-height (44 m).The non-dimensional time period on the horizontal axis corresponds to approximately one hour in actual time based on the wind velocity at the inflow boundary from the height H (=288 m) above the ground surface, U in = 7 m/s.(a) WTG-A.(b) WTG-B.

Figure 7 .
Figure 7.Comparison of the vertical profiles of the non-dimensional mean streamwise (x) wind velocity (symbol line), i.e., the mean streamwise wind velocity, U (= < u >), normalized by the streamwise velocity, U in , at H = 288 m above the ground surface at the inflow boundary.The vertical axis represents the height above the ground surface, z*, in full scale (m).The mean values were evaluated over a period of non-dimensional time of 100.The schematics in the figure illustrate the rotor diameter and the inflow velocity profile (solid line).(a) WTG-A; (b) WTG-B.

Figure 8 .
Figure 8.Comparison of the vertical profiles of the standard deviations of the three wind velocity components normalized by the streamwise (x) velocity, U in , at H = 288 m above the ground surface at the inflow boundary.The vertical axis represents the height above the ground surface, z*, in full scale (m).The standard deviations were evaluated over a period of non-dimensional time of 100.The schematics in the figure illustrate the rotor diameter.(a) WTG-A; (b) WTG-B.

Figure 9 .
Figure 9. Generation mechanism of turbulence (terrain-induced turbulence) contrast, at WTG-B, the maximum value of the standard deviation of the streamwise wind velocity component occurs near the hub-height of 44 m above the ground surface (Figure 8(b)).This height corresponds roughly to the height with large velocity shear as discussed for Figure 5(b) and Figure 7(b).
Fluid Dynamics and in the vicinity of WTGs will additionally be required.With the increasing size of WTGs, visualizations of vertical cross-sections of the flow field and evaluations of the vertical profiles of the mean wind velocity and turbulence intensity are expected to become increasingly important with time.
a turbulent boundary layer, in other words, this allows a simulation of turbulence in the incoming flow to a WF.The results from the present simulation set-up are discussed with a focus on the influence of the turbulence generated by the terrain near the WTGs on the WTGs.