Recently, the wind power industry has undergone rapid growth at an unprecedented rate across the world. This growth has been occurring because wind power generation has the best cost performance of all the renewable energies in terms of achieving a post-fossil fuel society and reducing CO2 emissions. There is no doubt that wind power is the leading renewable energy even in Japan. The authors are convinced that further dissemination of wind power generation will contribute on the global scale to “green innovation”, efforts to combat global warming.
In the field of wind power generation, as is the case in other fields, the use of CFD (Computational Fluid Dynamics) has increased rapidly and is being used for wind turbine deployment planning and estimation of the annual energy production of wind turbines. Given this background, validation testing of the prediction accuracy of CFD software used in the field of wind power generation has been advanced by various interested groups. One such validation testing was conducted as part of the Bolund Experiment  -  (see Figure 1). In this project, which was led by Risø DTU (Technical University of Denmark), the flow field passing over an isolated topographical feature protruding above the sea surface was targeted. During this project, wind condition observations (by cup and ultrasonic anemometers and LIDAR) and numerical simulations (CFD simulations) were performed; furthermore, based on these results, the prediction accuracy of CFD simulations has been analyzed. In these analyses, RANS (Reynolds Averaged Navier-Stokes) models (time-averaged models), which are steady, nonlinear wind analysis models, have received more favorable evaluations than LES (Large-Eddy Simulation) models (spatial-averaged models), which are unsteady, nonlinear wind analysis models. In contrast, in Japan, LES models as well as RANS models have received favorable evaluations. Accordingly, the present study re-examines the accuracy of airflow predictions by an LES model. For this purpose, the author conducted an independent wind tunnel experiment for the airflow over the topography studied in the Bolund Experiment. The present study compares the results from the LES-based numerical wind synopsis prediction technique RIAM-COMPACT, which was developed and is continuing to be developed by the author’s research group, and the airflow measurement results from the wind tunnel experiment.
2. Overview of the RIAM-COMPACT Natural Terrain Version Software
The core technology of RIAM-COMPACT (Research Institute for Applied
Figure 1. Bolund hill.
Mechanics, Kyushu University, COMputational Prediction of Airflow over Complex Terrain) is under continuous development at the Research Institute for Applied Mechanics, Kyushu University  -  . An exclusive license of the core technology has been granted by Kyushu TLO Co., Ltd. (Kyushu University TLO) to RIAM-COMPACT Co., Ltd. (http://www.riam-compact.com/), a venture corporation which was founded by the author and others and originated at Kyushu University in 2006. (A trademark, RIAM-COMPACT, and a utility model patent were granted to RIAM-COMPACT Co., Ltd. in the same year.) In the meantime, a software package has been developed based on the above-mentioned technique and is named the RIAM-COMPACT natural terrain version software. As of today, the RIAM-COMPACT software has been used by a large number of corporations and institutions including Eurus Energy Holdings Corporation and Electric Power Development Co., Ltd. (J-POWER) which have the largest share of the wind power generation industry in Japan.
Computation time had been an issue of concern for the RIAM-COMPACT software, which focuses on unsteady turbulence simulations. However, the present fluid simulation solver is compatible with multi-core CPUs (Central Processing Units) such as the Intel Core i7, which has drastically reduced the computation time, leaving no appreciable problems in terms of the practical use of the RIAM-COMPACT software. Furthermore, the RIAM-COMPACT software has successfully been made compatible with GPGPU (General Purpose computing on Graphics Processing Units). The concept of GPGPU is to widely apply the floating-point operation capacity of a GPU not only to graphics rendering but also to other numerical operations.
Subsequently, the numerical technique used for RIAM-COMPACT is described. Collocated grids in a general curvilinear coordinate system are used for the arrangement of variables in order to numerically predict local airflows over complex terrain with high accuracy while avoiding numerical instability. In these collocated grids, the velocity components and pressure are defined at the grid cell centers, and variables which result from the contravariant velocity components multiplied by the Jacobian are defined at the cell faces. As for the numerical method, the FDM (Finite-Difference Method) is adopted, and an LES model is used for the turbulence model. In LES, a spatial filter is applied to the flow field to separate eddies of various scales into GS (Grid Scale) components, which are larger than the computational grid cells, and SGS (Sub-Grid Scale) 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 relying on the use of a physically simplified model. On the other hand, 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 (Equation (1)) and a spatially filtered Navier-Stokes equation (Equation (2)) are used:
Supporting equations are given in Equations (3)-(8):
Because mean wind speeds of approximately 10 m/s or higher are considered in the present study, the effect of vertical thermal stratification (density stratification), which is generally present in the atmosphere, is neglected. The computational algorithm and the time marching method are based on a FS (Fractional-Step) method  and the Euler explicit method, respectively. The Poisson’s equation for pressure is solved by the SOR (Successive Over-Relaxation) method. For discretization of all the spatial terms in Equation (2) 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. An interpolation technique based on four-point differencing and four-point interpolation by Kajishima  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, α = 3.0 is commonly applied in the Kawamura-Kuwahara Scheme  . However, α = 0.5 is used in the present study to minimize the influence of numerical diffusion. For LES subgrid-scale modeling, the commonly used Smagorinsky model  is adopted. A wall-damping function is used with a model coefficient of 0.1.
3. Terrain Elevation Data and Simulation Set-Up
In this section, the terrain elevation data and simulation set-up are described. Figure 2 shows a contour map created from the terrain elevation data adopted in the present study. As indicated in this figure, DEM data with a spatial resolution of 1 m were created using and based on the terrain data provided by Risø DTU as the original data. The wind direction considered in the simulation is
Figure 2. Terrain elevation data used in the present study.
239˚, where the wind direction is defined such that 0˚ indicates northerly wind.
Figure 3 shows the evaluation points for the vertical profiles of airflow properties. In the present study, at four points (M1, M2, M3, and M4) on Line A which are defined for the study site of the Bolund Experiment, profiles of the mean wind velocity obtained from RIAM-COMPACT and the wind tunnel experiment are compared. As shown in the figure, the mean wind velocity profile acquired for point M0 on Line A from the wind tunnel experiment is given at the inflow boundary in the numerical simulation.
Figure 4 illustrates the computational domain, the computational grid, and the coordinate system used in the present study. The dimensions of the computational domain are 38.55 H, 38.55 H, and 15.45 H for the streamwise (x), spanwise (y), and vertical (z) directions, respectively, where H (=11 m, actual scale) is the maximum terrain elevation within the computational domain. The vertical dimension of the computational domain is set identical to that of the wind tunnel in order to avoid blockage effect. The numbers of computational grid points in the x-, y-, and z-directions are set to 213, 213, and 65, respectively, which results in a total of approximately three million grid points. As for the spatial resolution, the grid points are spaced at an even interval of Δx = Δy = 2 m (= 0.18 H) in the x-y cross section; the grid points are spaced at uneven intervals of Δz = 0.038 m to 7.64 m (0.003H to 0.7 H) in the z-direction so that the grid cells are distributed with higher density in the vicinity of the ground surface.
For the inflow conditions, as previously described in the discussion of Figure 3, the mean wind velocity profile obtained at point M0 on Line A in the wind tunnel experiment is given at the inflow boundary. The characteristic wind velocity scale Uref adopted in the present study is the value of the wind velocity at
Figure 3. Evaluation points for the vertical profiles of airflow properties.
Figure 4. Computational domain, computational grid, coordinate system, and other relevant information.
the inflow boundary at the height of the maximum terrain elevation, H (=0.07 m, wind tunnel scale). The Reynolds number Re (=UrefH/ν), which is based on Uref (=6.08 m/s) and H (=0.07 m), is set to 2.8 × 104 in order to match the value from the wind tunnel experiment.
4. Overview of the Wind Tunnel Experiment
Subsequently, an overview of the wind tunnel experiment performed in the present study will be given. Figure 5 shows a schematic view of the wind tunnel facility. This facility is the boundary layer wind tunnel of the Institute of Industrial Science, University of Tokyo. The dimensions of the test section are 16.5,
Figure 5. Boundary layer wind tunnel used in the present study (Test section length: 16.5 m, width: 2.2 m, height: 1.8 m).
2.2 and 1.8 m in length, width, and height, respectively. (See http://venus.iis.u-tokyo.ac.jp/en/introduction/facility/fudo/index.html)
Measurements of the airflow field were made using a split-fiber probe, which is capable of detecting reversed flow (Figure 6). The sampling intervals and the duration of the measurements were set to 1 kHz and 60 s, respectively.
Figure 7 shows the terrain model created for the present study. The scale of the terrain model was 1/150, and the maximum terrain elevation H (11 m) was approximately 0.07 m in the wind tunnel experiment (Figure 7(a)). Figure 7(b) is a close-up view of the terrain model. Figure 8 shows a view of airflow measurements conducted in the wind tunnel (viewed from downstream). Spires and roughness blocks were placed upstream of the model, and a boundary layer flow, in which the vertical profile of the mean wind speed followed a 1/7 power law distribution, was generated.
5. Results and Discussions
Figure 9 shows temporal changes of the non-dimensional instantaneous streamwise (x) velocity component u* (=u/Uref) at the measurement point as shown in the figure above. In Figure 9, the horizontal axis indicates the non-dimensional time t* (=t/(H/Uref)). An examination of Figure 9 reveals that a periodical flow phenomenon is generated behind the hill, that is, there exists the periodical vortex shedding. Attention is focused on the Figure 10 shows flow visualization of the instantaneous (non-dimensional time: t* = 141.40 - 146.20) wind fields obtained from the numerical simulation (RIAM-COMPACT LES turbulence model). Specifically, the figure shows the distribution of the non-dimensional instantaneous streamwise (x) velocity component on Line A (see Figure 3). Airflows around the hill include the unsteady vortex shedding in the wake region. That is, the flow separates in the vicinity of the downstream edge of the hill, and the separated shearing layer rolls up and forms isolated eddies.
Figure 6. Split-fiber probe, capable of detecting reversed flow.
Figure 7. Terrain model. (a) Overall view; (b) Enlarged view.
These small eddies then merge together. As a result, large vortices are formed, and are shed away to the downstream area of the model (arrow A in the figure).
Next, comparisons between the results from the numerical simulation (RIAM-COMPACT LES turbulence model) and those from the wind tunnel experiment are discussed.
Figure 11 shows visualization of the time-averaged (non-dimensional time: t* = 100 - 200) wind fields obtained from the numerical simulation (RIAM-COMPACT LES turbulence model). Specifically, the figure shows the distribution of the non-dimensional streamwise (x) velocity component. For comparison, the case in which the surface roughness was taken into consideration is also shown. The effect of the surface roughness was included by adding an external force term (a canopy model) to the governing equations (drag coefficient Cd = 10). This roughness model was applied to the air layer extending from the
Figure 8. View of airflow measurements conducted in the wind tunnel, viewed from downstream.
Figure 9. Results of the numerical simulations (RIAM-COMPACT LES turbulence model). Temporal changes of the non-dimensional instantaneous streamwise (x) velocity component at the measurement point as shown in the figure above (Non-dimensional time: t* = 100 - 200).
ground surface to z* = 0.0286 H (2 mm in the wind tunnel), where H is the maximum terrain elevation (11 m; 0.07 m in the wind-tunnel experiment), the characteristic length scale in the present study. This range for applying the roughness model was selected based on the results of the wind tunnel experiment. Examinations of the case with a smooth surface (Figure 11(a)) and the case with a rough surface (Figure 11(b)) reveal that vortex regions are formed downstream of the hill. Figure 11(b) clearly indicates elongation of the vortex region due to the effect of the surface roughness.
In Figure 12, the vertical profiles of the non-dimensional streamwise (x)
Figure 10. Results of the numerical simulations (RIAM-COMPACT LES turbulence model). Distribution of the non-dimensional instantaneous streamwise (x) velocity component. (Non-dimensional time: t* = 141.40 - 146.20). (a) Non-dimensional time: t* = 141.40; (b) Non-dimensional time: t* = 142.60; (c) Non-dimensional time: t* = 143.80; (d) Non-dimensional time: t* = 145.00; (e) Non-dimensional time: t* = 146.20.
Figure 11. Results of the numerical simulations (RIAM-COMPACT LES turbulence model). Distribution of the non-dimensional time-averaged streamwise (x) velocity component. (Non-dimensional time: t* = 100 - 200). (a) Case with no surface roughness (Smooth surface); (b) Case with surface roughness (Rough surface).
wind velocity at points M1 to M4 indicated in Figure 11 are compared. In this figure, lines, circles (○), and plus signs (+) indicate the results from the wind tunnel experiment (smooth surface), the LES numerical simulation (smooth surface), and the LES numerical simulation (rough surface), respectively. At all
Figure 12. Comparison between the numerical simulations (RIAM-COMPACT LES turbulence model) and the wind tunnel experiment. Plotted for the four points M1 to M4 shown in Figure 11. Vertical profiles of the non-dimensional streamwise (x) wind velocity component.
the examined points and heights, the results from the LES numerical simulations and wind tunnel experiment are in good agreement. Furthermore, at points M1 to M3, no significant differences due to the absence or presence of the surface roughness can be identified in the profiles of the mean wind velocity obtained from the LES numerical simulations. As described in the discussion of Figure 11, the effect of the surface roughness is evident in the size of the vortex region formed downstream of the hill. As a result of this effect, in the vicinity of the ground surface at point M4, there exist differences in the values of the mean wind velocity in the vertical profiles between the case with a smooth surface and the case with a rough surface.
In the present study, wind conditions were numerically predicted for the site of the Bolund Experiment using the RIAM-COMPACT natural terrain version software, which is based on an LES turbulence model (CFD). In addition, airflow measurements were made using a split-fiber probe in the boundary layer wind tunnel of the Institute of Industrial Science, University of Tokyo. The characteristics of the airflow at and in the vicinity of the site of the Bolund Experiment were clarified. The study also examined the prediction accuracy of the LES turbulence simulations (CFD). The values of the streamwise (x) wind velocity predicted by the CFD model were generally in good agreement with those from the wind tunnel experiment at all points and heights examined, demonstrating the validity of CFD based on LES turbulence modeling.
This work was supported by JSPS KAKENHI Grant Number 17H02053. Also, a wind tunnel experiment was supported by Dr. Makoto IIDA (Tokyo University). The author expresses appreciation to them.