Computational Fluid Dynamics ( CFD ) Approach to Predict an Actual Wind Speed over Complex Terrain

This paper proposes a procedure for predicting the actual wind speed for flow over complex terrain with CFD. It converts a time-series of wind speed data acquired from field observations into a time-series of actual scalar wind speed by using non-dimensional wind speed parameters which are determined beforehand with the use of CFD output. The accuracy and reproducibility of the prediction procedure were examined by simulating the flow with CFD with the use of high resolution (5 m) surface elevation data for the Noma Wind Park in Kagoshima Prefecture, Japan. The errors of the predicted average monthly wind speeds relative to the observed values were less than approximately 20%.


Introduction
A significant reduction in CO2 emissions to prevent global warming is currently an urgent issue.For this reason, the effective use of clean, environmentally-friendly wind energy has received attention.Since the power output of a wind turbine (WT) is proportional to the cube of the wind speed, it is important to accurately select a site with wind conditions that are well-suited for WT deployment with a high degree of precision.In Japan, most terrain is complex and few areas are flat, which is quite different from the terrain in Europe and the U.S.A.Thus, it is very important to take into account the effects of terrain on the wind flow, such as flow impingement, flow separation, flow reattachment, and reverse flow.Given this background, micro-siting software which can simulate wind over the complex terrain that is unique to Japan is currently being developed based on computational fluid dynamics (CFD).This development is being carried out with a vision that the software can also be applied for micro-siting of WTs over complex terrain in any other part of the world [1][2][3][4][5][6][7][8] .
RIAM-COMPACT is based on large eddy simulation (LES), which is a turbulence model, and simulates airflow in small domains of approximately a few hundreds of meters to a few kilometers.One of the major features of RIAM-COMPACT is that it can, with high accuracy, predict and visualize, in the form of animation, the effects of terrain on the wind flow such as flow separation, formation of reverse flow regions which results from flow separation, local flow acceleration, and reattachment of separated shear layers.
At a planned construction site for a wind power facility, it is a standard practice for an observation pole of 30 m or taller to be set up to collect time-series data of wind direction, speed, and other meteorological variables over a period of one year.If these observation data can be loaded into a CFD model and both the observation data and the CFD output can be combined efficiently, economic assessments of the wind power facility, including assessments of the annual energy production and capacity factor at an arbitrary location within the facility, can be made with high accuracy.
Accordingly, the present study proposes a procedure to predict the actual wind speed over complex terrain with the use of time-series data collected from field observations and CFD output.In this procedure, a time-series of wind speed data acquired from field observations is converted into a time-series of the actual scalar wind speed for the prediction point with the use of (non-dimensional) wind speed parameters which were determined beforehand with the use of CFD output.The validity of the proposed procedure is examined for Noma Wind Park, Kagoshima Prefecture, Japan, with the use of observation data collected with propeller vane anemometers mounted on the nacelles of WTs.

Overview of Noma Wind Park, Kagoshima Prefecture
Noma Wind Park is in Kasasa Town in the southwestern part of Kagoshima Prefecture, Japan (marked by the circle in Fig. 1).The terrain at and in the surrounding area of Noma Wind Park is typical complex terrain (Fig. 2).While the wind park is surrounded by the ocean, a steep cliff with a slope angle of more than 30 degrees is at the western end of the cape.The maximum elevation of the terrain surface is 143 m.In Noma Wind Park, ten WTs that are owned by Kyushu Electric Power Co., Inc. have been deployed, and verification tests have been conducted.The rated output of each WT is 300 kW, thus the total rated output from the ten WTs combined is 3000 kW.Tables 1 and 2 give an overview of the WTs deployed in Noma Wind Park.

Numerical Simulation Technique and Set-up
In order to numerically predict local wind flow over complex terrain with high accuracy while avoiding numerical instability, the present study employs RIAM-COMPACT.This software simulates wind flow with the use of a collocated grid in a general curvilinear coordinate system.In 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 simulation method, a finite-difference method (FDM) is adopted, and an 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.( 1)) and filtered Navier-Stokes equations (Eq.( 2)) are used.Because high wind conditions with mean wind speeds of 6 m/s or higher are considered in the present study, the effect of vertical thermal stratification, which is generally present in the atmosphere, is neglected.The effects of the surface roughness have been taken into consideration by reconstructing surface irregularities in high resolution (horizontal resolution: 5 m).
( ) For the computational algorithm, a fractional step (FS) method 19) is used, and a time marching method based on the Euler explicit method is adopted.The Poisson's equation for pressure is solved by a successive over-relaxation (SOR) method.For discrimination of all the spatial terms in the governing equations except for the convective term in Eq. ( 2), a second-order central difference scheme is applied.For the convective term, a third-order upwind difference scheme is used.For the fourth-order central difference term that constitutes the discretized convective term, the interpolation technique of Kajishima 20) , which is based on a four-point difference scheme and a four-point interpolation scheme, is used.For the weighting of the numerical diffusion term in the discretized convective term, α = 3.0 is commonly applied in the Kawamura-Kuwahara scheme 21) .However, α = 0.5 is used in the present study to minimize the influence of numerical diffusion.For the LES subgrid-scale modeling, the standard Smagorinsky model 22) (Eqs.( 3) -( 8)) is adopted with a model coefficient of 0.1 in conjunction with a wall-damping function.

Procedure for Predicting Actual Wind Speed
The procedure for predicting the actual scalar wind speed time-series data for an arbitrary point over complex terrain is described below.In this procedure, time-series data that were acquired from field observations and CFD output (time-averaged wind speed) are utilized.
Step 1.As illustrated in Fig. 6, wind simulations are performed by RIAM-COMPACT for 16 individual wind directions for a prediction point.After the flow field has sufficiently developed in a simulation run, the simulation is continued, and the time-averaged flow is calculated.The computational domain is a square with 5 km sides with WT #4 in the center of the domain (Fig. 6).The top of the computational domain is set to 700 m.Terrain features are reconstructed with the use of high-resolution (5 m spatial resolution) surface elevation data.The numbers of computational grid points are 51 (streamwise (x)) × 51 (spanwise (y)) × 41 (vertical (z)).Variable grid spacing is used in the x-and y-directions so that the grid density increases toward the center of the domain.Likewise, variable grid spacing is used in the z-direction so that grid density is highest near the ground surface.
The vertical height of the first layer of grid cells, that is, the minimum vertical grid height, is approximately 2.5 m.A vertical profile of wind speed based on the 1/7 power law is given for the inflow boundary.The free-slip condition is set for the side and upper boundaries.The convective outflow condition is applied to the outflow boundary.On the ground, the no-slip condition is imposed.The non-dimensional parameter, Re, in Eq. ( 2), defined as Uin h/ν is the Reynolds number and is set to Peer-reviewed version available at Energies 2018, 11, 1694; doi:10.3390/en11071694 Step 2. For each inflow direction (Fig. 6), the ratio of the simulated hub-height wind speed at a prediction point to the simulated wind speed at the height of an observation pole, that is, at the reference point, is computed, yielding a set of 16 wind speed ratios for the prediction point.The value of the wind speed at the reference point is normalized as "1."In this step, the wind direction set at the inflow boundary is assumed to match with that at the reference point.In the present study, points located at WT #4 and WT #6 at the hub-height are used as the reference points.(Hereafter, in the context of reference points, these reference points will be referred to as WT #4 and WT #6.)Thus, the following wind speed ratios are evaluated: WT #1 / WT #4, WT #2 / WT #4, WT #3 / WT #4, WT #5 / WT #4, …, WT #10 / WT #4; WT #1 / WT #6, WT #2 / WT #6, WT #3 / WT #6, WT #4 / WT #6, WT #5 / WT #6, WT #7 / WT #6, …, WT #10 / WT #6.To compute these wind speed ratios, the scalar wind speed, VEL, is used in both the denominator and numerator and is computed from the time-averaged wind velocity components in the x-and y-directions as ( ) Step 3. The observed time series of the scalar wind speed at the reference point is converted to a time-series of scalar wind speed at a prediction point as follows.According to the wind direction at the reference point at a given instance in the observed time series, the scalar wind speed at the reference point is multiplied by the wind speed ratio which was determined beforehand for the prediction point for that wind direction in Step 2. By applying this operation to an observed time series of scalar wind speeds that extends over a period of one year or any period of interest and processing the outcome of the operation statistically, the monthly and annual average wind speed can be obtained for an arbitrary prediction point.predicted time-series of scalar wind speed for the hub-height of WT #1 and WT #10, however, WT #6 was used as the reference point for the prediction of the scalar wind speed.Scatter plots for the data in Fig. 11 are shown in Fig. 12.In both the case of using WT #4 and the case of using WT #6 as the reference points, the proposed procedure for predicting the actual wind speed yielded predicted values which agrees well with the observed values.

Results and Discussions
Table 3 summarizes comparisons of observed and predicted values of wind speed for the cases of adopting WT #4 and WT #6 as the reference points.Table 4 gives the horizontal separation distance between a given WT and WT #4 or WT #6, both of which were used as reference points.Similarly, Table 5 shows the elevation difference between the hub-height of a given WT and reference point WT #4 or WT #6.Examinations of Tables 3 to 5 reveal that the relative error of the wind speed at a WT predicted by the proposed actual wind speed prediction procedure is less than approximately 20% if the horizontal separation distance and elevation difference between a prediction point (the hub-height of a WT) and the reference point are approximately 800 m or less and 50 m or less, respectively.As expected, the observed wind speed data that are used in the present study were affected by the blade rotation since they were collected from propeller vane anemometers that were mounted on the nacelles of the WTs.However, this effect is equally present in the wind speed data that were collected at all the WTs.In regard to the previously discussed range of application for the proposed actual wind speed prediction procedure, the subsequent analysis examines, based on the observation data, the characteristics of airflow over complex terrain, which serve as supporting evidence for the range of application for the procedure.In this analysis, an angle of 0° denotes wind which blows from north to south, and wind which flows from wind directions that are 348.75°or larger and 11.25°or smaller are defined as northerly wind.Fig. 14 shows the time-series of the wind speed ratio only for times of northerly wind.Fig. 14 also shows the average value of the observed wind speed ratios for northerly wind and that from the CFD model, revealing that these two average values agree closely.Fig. 15 shows the frequency distribution of the wind speed ratios plotted in Fig. 14.The frequency is distributed nearly symmetrically with the average wind speed ratio in the center.Thus, it can be said that time-series data of scalar wind speed at an arbitrary point can be obtained with the use of observed time-series of wind speed from the reference point if, with the use of CFD, the effects of terrain irregularities on the wind flow can be simulated accurately and the ratio of the wind speed at the prediction point to the wind speed at the reference point can be estimated with high accuracy.Peer-reviewed version available at Energies 2018, 11, 1694; doi:10.3390/en11071694

Conclusion
In the present study, a procedure for predicting the actual wind speed over complex terrain was proposed.This procedure utilizes both time series of observation data and CFD output obtained from the unsteady, nonlinear wind simulator RIAM-COMPACT, in which a collocated grid in a general curvilinear coordinate system is adopted.In this procedure, time-series of scalar wind speed obtained from field observations are converted to the actual scalar wind speed with the use of (non-dimensional) wind speed parameters that were obtained beforehand with the use of the CFD.The proposed prediction procedure was validated with the use of observation data acquired with propeller vane anemometers mounted on the nacelles of wind turbines in Noma Wind Park, Kagoshima Prefecture, Japan.As for the result, the following was found for Noma Wind Park, which was investigated in the present study.If the horizontal separation distance between a prediction point and the reference point is approximately 800 m or less and the elevation difference between the two locations is 50 m or less, the relative error of the wind speed predicted by the proposed wind speed prediction procedure is less than approximately 20%.
However, additional validation is necessary in regard to how the reference location is selected and the types of observation instruments used at the reference point.In addition, validation also needs to be carried out at other sites.These validations will be topics of future research.Finally, if the proposed procedure is further expanded so that it utilizes output from a meteorological model instead of observation data, a so-called real-time simulation, which can predict wind flows a few hours to a few days into the future, will be made possible.

Fig. 6
Fig. 6 Conceptual outline of wind simulations for 16 wind directions for Noma Wind Park

Fig. 8 Fig. 8 1 (b) WT # 10 Fig. 9 10 Fig. 10 4 (Fig. 11 1 (b) WT # 10 Fig. 12 6 Fig. 9
Fig.8illustrates instantaneous views of flows that were visualized with the passive particle tracking method.The inflow wind direction for this case is northerly.This figure indicates that, due to the influence of the terrain irregularities, wind speed and direction change in complex ways as the wind flows over Cape Noma.In the procedure for predicting the actual wind speed which is proposed in the present study, simulations are run until fully-developed flows, as the one illustrated in Fig.8, are achieved, and the simulations are continued for a period of time.Then, the results of the simulations are used to obtain the time-averaged flow field.

Fig. 13
Fig. 13 Comparison of hub-height wind speed at WT #1 and WT #10, observation data Fig. 13 shows a comparison of the time-series of the observed wind speed from WT #1 and WT #10 and also a scatter plot of the observed wind speed from the two WTs for February 2003.The horizontal separation distance between the two WTs is as much as approximately 1200 m, and the elevation difference between the two WTs is 9 m (see Fig. 2).In particular, attention should be given to the horizontal Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 24 May 2018 doi:10.20944/preprints201805.0328.v1Peer-reviewed version available at Energies 2018, 11, 1694; doi:10.3390/en11071694separation distance.Even though the horizontal separation distance between WT #1 and WT #10 is more than 1000 m, the temporal changes of the wind speed over the course of one month at WT #1 and WT #10 are quite similar.The correlation coefficient between the two data sets plotted in Fig. 13(b) is 0.84.Of the 40,047 values of the ratio of the wind speed at WT #1 (the prediction point) to the wind speed at WT #4 (the reference point) for February 2003, data values only from the time of northerly wind are extracted in order to study the statistical properties of the wind speed ratio.The number of extracted data values is 9,501, which corresponds to approximately 24% of the entire set of data values.

Fig. 14
Fig. 14 Observed time-series of the ratio of the hub-height wind speed at WT #1 to that at WT #4, only for times of northerly wind, February 2003

Table 1
Specifications of WTs at Noma Wind Park

www.preprints.org) | NOT PEER-REVIEWED | Posted: 24 May 2018 doi:10.20944/preprints201805.0328.v1
View of Noma Wind Park, the numbers in the figures are those of the WTs

Table 3
Comparisons of observed and predicted values of wind speed (a) Reference point: WT #4

Table 4
Horizontal separation distance between a given WT and the reference point

Table 5
Elevation difference between the hub-height of a given WT (a prediction point) and the reference point