Comparison of RANS and LES in the Prediction of Airflow Field over Steep Complex Terrain

The present study compared the prediction accuracy of the three CFD software packages for simulating airflow around a three-dimensional, isolated hill with a steep slope: 1) WindSim (turbulence model: RNG k-ε RANS), 2) Meteodyn WT (turbulence model: k-L RANS), which are the leading commercially available CFD software packages in the wind power industry and 3) RIAM-COMPACT (turbulence model: standard Smagorinsky LES), which has been developed by the lead author of the present paper. Distinct differences in the airflow patterns were identified in the vicinity of the isolated hill (especially downstream of the hill) between the RANS results and the LES results. No reverse flow region (vortex region) characterized by negative wind velocities was identified downstream of the isolated hill in the result from the simulation with WindSim (RNG k-ε RANS) and Meteodyn WT (k-L RANS). In the case of the simulation with RIAM-COMPACT natural terrain version (standard Smagorinsky LES), a reverse flow region (vortex region) characterized by negative wind velocities clearly forms. Next, an example of wind risk (terrain-induced turbulence) diagnostics was presented for a large-scale wind farm in China. The vertical profiles of the streamwise (x) wind velocity do not follow the so-called power law wind profile; a large velocity deficit can be seen between the hub center and the lower end of the swept area in the case of the LES calculation (RIAM-COMPACT).


Introduction
We have developed an unsteady and non-linear wind synopsis simulator called RIAM-COMPACT (Research Institute for Applied Mechanics, Kyushu Univer-sity, 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]- [13].In RIAM-COMPACT, a large-eddy simulation (LES) has been adopted for turbulence 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.Efforts have been made to promote RIAM-COMPACT, mainly in the wind power industry (e.g., private wind power providers, local governments, and wind turbine manufacturers) in Japan.Computation time had been an issue of concern for the RIAM-COMPACT software, which focuses on unsteady turbulence simulations (LES).The present fluid simulation solver is compatible with multi-core CPUs such as the Intel Core i9 and also with GPGPU, which has drastically reduced the computation time, leaving no appreciable problems in terms of the practical use of the RIAM-COMPACT software.
On another front, commercially available CFD software such as STAR-CCM+ [14] and ANSYS (CFD, Fluent, CFX) [15] has developed mainly as an engineering tool primarily in the automobile and aviation industries until the present time.Recently, some of the above-mentioned general purpose thermal fluid analysis software has started being adopted in the wind power industry.In the previous study [16] [17] [18], the simulation results obtained from the RIAM-COMPACT software were compared to those from STAR-CCM+, one of the leading commercially available CFD software packages.The results of the comparison are discussed.In addition, open-source CFD software packages are more widely used than in the past.One of the most widely used software packages is OpenFOAM (OpenField Operation And Manipulation) [19].Open-FOAM is an open-source CFD toolbox which has been released and distributed under the GNU GPL (General Public License) [20] by the OpenFOAM Foundation, a non-profit organization.In the previous study [21], the simulation results obtained from the RIAM-COMPACT software were also compared to those from OpenFOAM, and the results of the comparison are discussed.
The wind power industry has on its own developed and distributed CFD software designed for selecting sites appropriate for the installation of wind turbine generators.One such leading software package is Meteodyn WT [22], which has been developed by Meteodyn in France.Meteodyn WT is a CFD software package which incorporates a RANS turbulence equation with a one-equation closure scheme (k-L turbulence model; here, k and L refer to turbulence energy and the turbulence length scale, respectively).On October 12, 2017, "WT6.0", the latest version of the software, was released.Another one of the most widely used software packages is WindSim [23] by Norway-based WindSim AS.WindSim is a CFD software package which uses a RANS turbulence model and has been designed specifically for wind resource assessment.In December 2016, the latest version of the company's software package, WindSim 8.0, was released.These two CFD software packages are specialized for wind power resource assessment as well as the RIAM-COMPACT CFD software package.
In the present study, numerical simulations are performed for airflow over and around a three-dimensional, isolated hill with a steep slope angle using the three CFD software packages (WindSim, Meteodyn WT and RIAM-COMPACT).The results of the comparison are discussed.Next, the numerical simulations for airflow over a large-scale wind farm in China [21] are performed with RIAM-COMPACT, which is based on an LES turbulence model, and Meteodyn WT, which is based on a RANS turbulence model.The numerical wind simulations in the present study are conducted for high Reynolds number airflow over and around a three-dimensional, isolated hill with a steep slope angle and a large-scale wind farm in China.Table 1 shows the simulation set-ups adopted in the two software packages which are used in    1 and Figure 2 illustrate the full computational grid used for WindSim and an enlarged view of the grid in the vicinity of the isolated hill, respectively.Figure 3 shows the inflow profile used for all of the simulations in the present study.Figure 4 shows the characteristic wind velocity and length scales which are employed for the simulations with RIAM-COMPACT.

Overview of the
In RIAM-COMPACT, a collocated grid in a general curvilinear coordinate system is used in order to numerically predict local wind flow over complex terrain with high accuracy while avoiding numerical instability.For the numerical simulation method, a FDM is adopted, and an LES model is used for the turbulence model.For the computational algorithm, a method similar to a FS method [24] 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 SOR method.For discretization of all the spatial terms in the governing equations except for the convective term in the Navier-Stokes equation, a second-order central difference scheme is applied.For the convective term, a third-order upwind difference scheme is used.An interpolation technique based on four-point differencing and four-point interpolation by Kajishima [25] 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 [26].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 [27] is adopted with a model coefficient of 0.1 in conjunction with a wall-damping function.For further details of the numerical simulation techniques, refer to Uchida [1]- [13].
Regarding the boundary conditions adopted for the simulations with RIAM-COMPACT, the same inflow profile as used for the simulations with WindSim (Figure 3) is given at the inflow boundary.At the side and upper boundaries, free-slip conditions are applied, and convective outflow conditions Open Journal of Fluid Dynamics are applied at the outflow boundary.On the ground surface, a non-slip boundary condition is imposed.For the simulation at Re (=U in h/ν) = 10 7 , the number of grid points is changed to 101 in the vertical direction, and the minimum vertical grid spacing in is set to Δz min /h = 4 × 10 −7 according to the equation below (see Table 1): In contrast to RIAM-COMPACT, WindSim uses RANS models.In the present study, the RNG k-ε RANS model is selected for the simulations.Refer to [23] for the numerical simulation methods used in WindSim and other details about this software.Figure 9. Method adopted in Meteodyn WT for setting the inflow profile and the inflow profile generated for the present study.

Comparison of the Simulation
Since simulations for a flow with Re (=U in h/ν) = 10 7 were not feasible with the RIAM-COMPACT natural terrain version software because of the time step, a numerical wind simulation is performed at Re = 10 6 , which is an order of magnitude smaller than the flow simulated with Meteodyn WT.For this simulation, the number of grid points in the vertical direction is set to 101 (37 for the simulation with Meteodyn WT), and the minimum vertical spacing is set to Δz min /h = 10 −4 based on the Equation (1) (Δz min /h = 5.0 × 10 −3 for the simulation with Meteodyn WT, see Table 2).At the inflow boundary, an inflow profile which is almost identical to the inflow profile used for the simulation with Meteodyn (Figure 9) is used.Free-slip conditions are applied at the side and upper boundaries, and convective outflow conditions are applied at the outflow boundary.At the surfaces of the ground and the isolated hill, non-slip conditions are imposed.The time step is set to Δt = 10 −5 h/U in (refer to Table 2).The results from the simulations are compared.

Comparison of the Simulation Results from the Two CFD Software Packages (RIAM-COMPACT and Meteodyn WT) in the Case of a Three-Dimensional, Isolated Hill with a Steep Slope Angle
The vibration problem of turbine T12 was investigated by the operator Open Journal of Fluid Dynamics YUDWPC and a report was issued in April 2014 [28].Stated in the report was that high vibration data was recorded only when the wind was blowing from the southwest.Wind direction on the ground level was observed to be in the reverse direction from that recorded by the nacelle anemometer.Analysis of the vibration data indicates the vibration is in the vertical direction.This suggests the vibration is associated with abnormal vertical wind shear across the wind turbine rotor.As shown in Figure 15, a figure extracted from the report, it was deduced that the presence of the small hill located about 150 m upstream from turbine T12 was causing the onset of turbulence and reverse flow which led to the vibration recorded.
For LES simulation, the RIAM-COMPACT natural terrain version software package was employed.The software uses a standard Smagorinsky turbulence model.For the simulation, SRTM 90 m data was used for elevation data.Wind direction is set to true north at 247 degrees and the computational domain constructed is shown in Figure 16 with the following details:  To increase calculation accuracy, the mesh is concentrated around the turbine positions in both the x and y direction as shown in Figure 16.No roughness consideration is given in the present LES simulation.Atmospheric stability is set to neutral stability.After the calculation has been stabilized, numerical results in the calculation domain are output for a real time of ten minutes with an interval of one second.In this study, the commercial software Meteodyn WT (turbulence model: k-L RANS) was employed and its results were compared with the results calculated by the RIAM-COMAPCT.The calculation parameters are shown in Table 3.
The calculation domain is a radius of 10 km for the x-y direction with turbine T12 as center; z direction has a maximum of 200 m.Wind Direction is set at 247 degrees with minimum vertical and horizontal resolution set to 5 m and 2 m respectively.Atmospheric stability is set to neutral.The calculation was completed smoothly with computation convergence recoded at 99.3%.20 and the resulting wind shear profile is compared with the shear profile (average values) predicted by RIAM-COMPACT.It can be seen from Figure 20 that the shapes of the two profiles are similar from 50 m upwards but distinctively different below 50 m.Meteodyn WT does not seem to predict any flow separation and reverse flow region and therefore there is no significant wind speed reduction between 25 m and 50 m, and also no negative wind speed values below 25 m as predicted by RIAM-COMPACT.Numerical comparison results are shown in Table 5.
Referring to Table 5, across the wind turbine rotor face, RIAM-COMPACT predicted a large wind speed difference with a shear exponent exceeding the IEC standard value of 0.2 by a large margin.In sharp contrast, Meteodyn WT predicted a small wind speed difference with a shear exponent of 0.025 which is significantly below the IEC standard.

Summary
Simulations were performed for airflow around a three-dimensional, isolated hill with a steep slope angle in order to compare the flow pattern simulated in the vicinity of the hill by three software packages.For the simulations, three software packages were used: 1) WindSim (turbulence model: RNG k-ε RANS), 2) Meteodyn WT (turbulence model: k-L RANS), which are the leading commercially available CFD software packages in the wind power industry and 3) RIAM-COMPACT (turbulence model: standard Smagorinsky LES).Comparisons of the simulated results revealed a distinct difference in the simulated flow patterns in the vicinity of the isolated hill (especially downstream of the hill).No reverse flow region (vortex region) characterized by negative wind velocities was identified downstream of the isolated hill in the result from the simulation with WindSim (RNG k-ε RANS) and Meteodyn WT (k-L RANS).In the case  Thus, the flow pattern which forms in the vicinity of the isolated hill varies significantly according to the velocity boundary conditions applied for the surfaces of the ground and the isolated hill.
Software Packages (RIAM-COMPACT and WindSim) and Numerical Simulation Set-Up in the Case of a Three-Dimensional, Isolated Hill with a Steep Slope Angle

Figure 3 .
Figure 3. Inflow wind velocity profile used in the present study.

Figure 4 .
Figure 4. Characteristic wind velocity and length scales used in the simulation with RIAM-COMPACT.

Figure 5 .
Figure 5 shows the ensemble-averaged flow fields from the simulations with WindSim (turbulence model: RNG k-ε RANS).In neither of these simulations

Figure 6 Figure 6 .
Figure 6 shows instantaneous flow fields from the simulations with RIAM-COMPACT (turbulence model: standard Smagorinsky LES) (Re = 5 × 10 4 and 1 × 10 7 ).An examination of these simulation results reveals the clear presence of a reverse flow region (vortex region), in which the values of the streamwise wind velocity are negative, downstream of the isolated hill.

Figures 10 -
Figures 10-12 show results from the simulation with the Meteodyn WT software package (turbulence model: k-L RANS).These results (for a flow at Re = 10 7 ) indicate that a reverse flow region (vortex region) characterized by negative values of wind velocity does not form downstream of the isolated hill, and a pattern resembling potential flow is present there.Figure 13 shows the results from the simulation with the RIAM-COMPACT natural terrain version software package (turbulence model: the standard Smagorinsky LES).Examinations of the results reveal that a reverse flow region (vortex region) characterized by negative values of wind velocities clearly exists downstream of the isolated hill in the simulated flow at Re (=U in h/ν) = 10 6 .

Figure 12 .
Figure 12.Velocity vectors at the center of the span (y = 0) in the vicinity of the isolated hill, Meteodyn WT, k-L RANS model, Re = 10 7 .

Figure 13 . 6 .
Figure 13.Streamwise (x) wind velocity distribution at the center of the span (y = 0) in the vicinity of the isolated hill, RIAM-COMPACT, standard Smagorinsky LES, Re = 10 6 .(a) Instantaneous flow field; (b) Time-averaged flow field.

Figure 15 .
Figure 15.Deduction made on the airflow upstream and in the vicinity of turbine T12.

Figure 16 .
Figure 16.Computational domain and grid used for the simulation with RIAM-COMPACT.

Figure 17 showsFigure 17 .
Figure 17 shows an instantaneous vector plot across turbine T12.This picture clearly shows flow separation occurred at the small hill located 140 m upstream from the turbine, and the onset of the formation of the recirculating vortex behind the hill.The turbulent flow extends downstream forming a reverse flow region characterized by negative values of wind speed covering the lower part of the wind turbine rotor.The simulation results also indicate that the wind flow is relatively undisturbed above hub height level.The U component wind speed time series during the ten minute simulation at rotor top (106.3m), hub center (65 m), rotor bottom (23.7 m) and surface level (10 m) positions are plotted in Figure 18.Referring to Figure 18, it is obvious that the wind speed at surface (10 m) and rotor bottom is significantly lower and showing more fluctuations than the wind speed at the hub and top part of the rotor.Wind speed varies between 15.0 to 20.0 m/s at hub height and rotor top whereas for rotor bottom wind speed fluctuates between negative 6.2 m/s to 2.0 m/s.Negative wind speed indicates the

Figure 18 .
Figure 18.Time series of U component wind speed at rotor top, hub, rotor bottom and ground surface level at turbine T12.

Figure 19 .
Figure 19.Vertical shear profile predicted by RIAM-COMAPCT, average, minimum and maximum of U component wind speed variation with height.

Convergence 99. 3 %
Meteodyn WT's calculation output includes the speed-up factor from height 20 m to 200 m at an interval of 20 m at turbine T12.These values are shown in Table 4.The speed-up factor is the wind speed ratio at the given height referencing the wind speed at height 10 m.The speed-up factor therefore resembles the vertical shear profile.Assuming a wind speed of 9.5 m/s at 10 m height, wind

Figure 20 .
Figure 20.Comparison of vertical shear profile between the wind flows simulated by RIAM-COMPACT and Meteodyn WT.

Figure 22 .
Figure 22.Streamwise (x) wind velocity distribution in the vicinity of the isolated hill at the center of the span (y = 0), RIAM-COMPACT, standard Smagorinsky LES, Re = 10 4 .(a) Case 1: Simulation result from the case in which non-zero wind velocities are appliedas a Dirichlet boundary condition; (b) Case 2: Simulation result from the case in which all three wind velocity components are set to zero as a Dirichlet boundary condition.

Table 3 .
Numerical simulation methods, parameters, and settings for Meteodyn WT.

Table 4 .
Speed-up factor at turbine T12 and calculated wind speed.

Table 4 .
The wind speed figures in Table4are plotted in Figure