Reports of Research Institute for Applied Mechanics, Kyushu University No.149 (91–98) September 2015

Analysis of the Airflow Field around a Steep, Three-dimensional
Isolated Hill with Commercially Available CFD Software
Takanori UCHIDA*1, Fumihito WATANABE*2 and Shin MIKAMI*2
E-mail of corresponding author: takanori@riam.kyushu-u.ac.jp(Received August 31, 2015)

Abstract
The present study compared the prediction accuracy of two CFD software packages for simulating a flow
around a three-dimensional, isolated hill with a steep slope: 1) Meteodyn WT (turbulence model: k-L RANS), which
is one of the leading commercially available CFD software packages in the wind power industry and 2)
RIAM-COMPACT® (turbulence model: the standard Smagorinsky LES), which has been developed by the leadauthor of the present paper. Although the Reynolds number settings differed between the simulation with Meteodyn
WT and that with RIAM-COMPACT®, distinct differences in the flow patterns were identified in the vicinity of theisolated hill (especially downstream of the isolated hill) between the flows simulated by the two CFD software
packages.
Key words : Meteodyn WT, RIAM-COMPACT® , Isolated hill

1. Introduction
The lead author has developed the numerical winddiagnosis technique named RIAM- COMPACT®1),
which is based on an LES turbulence model. An
exclusive license for RIAM-COMPACT® has been
granted by Kyushu TLO Co., Ltd. to RIAM-COMPACTCo., Ltd. (http://www.riam-compact.com/), a venturecorporation which originated at Kyushu University.
Efforts have been made to promote this software mainlyin the wind power industry (e.g., private wind power
providers, local governments, and wind turbinemanufacturers) in Japan.
At the same time, other original CFD software
packages have been developed outside Japan and havebeen distributed in various countries. One such leading
software package is Meteodyn WT2), which has beendevelopedbyMeteodyn
inFrance(http://meteodyn.com/). Meteodyn WT is a CFD software
package which incorporates a RANS turbulence modeland is specialized for wind power resource assessment.
On January 20, 2014, "WT4.5", the latest version of thesoftware, was released.
In the present paper, numerical simulations areconducted for airflow around a three-dimensional,*1

Research Institute for Applied Mechanics,
Kyushu University*2Hitachi Power Solutions Co., Ltd

isolated hill with a steep slope angle using RIAMCOMPACT®, which is based on an LES turbulencemodel, and Meteodyn WT, which is based on a RANS
turbulence model. The computational results from thesesimulations are compared.

2. Overview of the two software packages
and the computational parameters
For the present study, numerical simulations areconducted for high Reynolds number flow around athree-dimensional, isolated hill with a steep slope angle
using RIAM-COMPACT®, which is based on an LES
turbulence model, and Meteodyn WT, which is based on
a RANS turbulence model (see Table 1). Fig. 1 shows
the computational domain and grid used for thesimulation with Meteodyn WT. An enlarged view of thecomputational grid in the vicinity of the isolated hillfrom the simulation with Meteodyn WT is shown in Fig.
2. Fig. 3 shows the method used to set the inflow
profile in Meteodyn WT and the inflow profile
generated for the present study. Fig. 4 illustrates thecharacteristic scales used in RIAM-COMPACT®; forthe present study, h is the height of the isolated hill, Uin
is the wind velocity at the inflow boundary at height h,and ν is the coefficient of dynamic viscosity.
The computational techniques and parameters used
for the simulations with the two software packages are

92

Uchida et al. : Analysis of the Airflow Field with Commercially Available CFD Software

CFD model

RIAM-COMPACT®

Turbulence model

METEODYN WTThe standard Smagorinsky

k-L RANS model

LES model(A single equation model)

Atmospheric stratification

Neutral(Atmospheric stability)
Coriolis force

Not included
Roughness length: 0.05

Surface roughness

Not included (Smooth surface)(For the ground surface
not on the isolated hill: 0.001)

Surface boundary condition
Shape function of the isolated hill

No-slip boundary condition(The three components of the wind velocity are set to zero.)
0.5h×{1+cos(πr/a)}, where r = (x2+y2)1/2, a = 2hHeight of the isolated hill h100 (m)Reynolds number Re (= Uinh/ν)

10107

Time step ∆t

10-5h/Uin (s)-Dimension of the computational
domainNumber of computational grid points

619h (i) × 18h (j) × 8h (k)436 (i) × 325 (j) × 101 (k)436 (i) × 325 (j) × 37 (k)(Approx. 14.3 million)(Approx. 5.2 million)

Streamwise (x) grid spacing (⊿x)(0.035 - 0.5) h

Spanwise (y) grid spacing (⊿y)
Vertical (z) grid spacing (⊿z)

(0.0001 - 0.6) h(0.005 - 1.2) hTable 1 Comparison of computational techniques and parameters
used for the simulations with RIAM-COMPACT® and Meteodyn WTApprox. 8hApprox. 9h
Approx. 13h
j = the 163rd grid point

Approx. 9h
Approx. 6h

Fig.1

i = the 113th grid point

Computational domain and grid used for the simulation with Meteodyn WT, Re = 107

Reports of Research Institute for Applied Mechanics, Kyushu University No.149 September 2015

Fig.2

93

Enlarged view of the computational grid used for the simulation with Meteodyn WT, Re=107

In Meteodyn WT, the inflow wind is not a uniform flow. The inflow windvelocity changes with height in order to represent actual atmosphericflow. Specifically, the inflow wind can be classified into the following threelayers.
• Surface layer: 0 < z < hs→The profile is determined by the logarithmic law.
• The Ekman layer: hs < z < hbl→Transition from the upper-surface-layer wind to the geostrophic wind inthe free atmosphere• Free atmosphere: hbl < z < Zskv
→The geostrophic wind is applied throughout this layer including as thetop boundary condition.

Normalized

The inflow profile is calculated based on the following reference values:•Vref = 10.0 m/s: the reference wind velocity•Href = 10.0 m: the reference height•Zoref = 0.05 m: the reference roughness length
In the surface layer, the inflow profile is calculated based on theatmospheric stability, surface roughness data from the vicinity of theinflow boundary, and an assumption that turbulence kinetic energy isconstant across the layer.

Isolated hill

In the Ekman layer, turbulence kinetic energy is assumed to decreaselinearly with height.

Fig.3

Method adopted in Meteodyn WT for setting the inflow profile andthe inflow profile generated for the present study

Uchida et al. : Analysis of the Airflow Field with Commercially Available CFD Software

listed in Table 1. The present study uses theRIAM-COMPACT® natural terrain version software
package, which is based on a collocated grid in a generalcurvilinear coordinate system, in order to numerically
predict local wind flow over complex terrain with highaccuracy while avoiding numerical instability. In thiscollocated 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 theLES model, a spatial filter is applied to the flow field toseparate 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 effectof small-scale eddies, i.e., the SGS components, ismodeled according to a physics-based analysis of theSGS stress.
For the governing equations of the flow, a filteredcontinuity equation for incompressible fluid and afiltered Navier-Stokes equation are used. Because highwind conditions with mean wind speeds of 6 m/s orhigher are considered in the present study, the effect ofvertical thermal stratification, which is generally present
in the atmosphere, is neglected (neutral atmosphere).
The effect of the surface roughness is not taken into
consideration, and the surfaces of the ground and theisolated hill are treated as smooth surfaces.
For the computational algorithm, a fractional step(FS) method is used, and a time marching method basedon the Euler explicit method is adopted. The Poisson’sequation for pressure is solved by the successive overrelaxation (SOR) method. For discretization of all thespatial terms except for the convective term in the
Navier-Stokes equation, a second-order central
difference scheme is applied. For the convective term, athird-order upwind difference scheme is applied. Aninterpolation technique based on four-point differencing
and four-point interpolation by Kajishima is used for thefourth-order central differencing that appears in thediscretized form of the convective term. For theweighting of the numerical diffusion term in theconvective term discretized by third-order upwinddifferencing, α = 3.0 is commonly applied in the

Kawamura-Kuwahara scheme. However, α = 0.5 is usedin the present study to minimize the influence ofnumerical diffusion. For LES subgrid-scale modeling,
the standard Smagorinsky model is adopted with amodel coefficient of 0.1 in conjunction with awall-damping function. Refer to Uchida1) for details of
the computational techniques.
Since simulations for a flow with Re (= Uinh/ν) = 107
were not feasible with the RIAM-COMPACT® natural
terrain version software because of the time step, anumerical wind simulation is performed at Re = 106,
which is an order of magnitude smaller than the flow
simulated with Meteodyn WT. For this simulation, thenumber of grid points in the vertical direction is set to101 (37 for the simulation with Meteodyn WT), and theminimum vertical spacing is set to ⊿Zmin/h = 10-4 basedon the equation below (⊿Zmin/h = 5.0×10-3 for thesimulation with Meteodyn WT, see Table 1):

z min / h 0.1
ReAt the inflow boundary, an inflow profile which isalmost identical to the inflow profile used for thesimulation with Meteodyn (Fig. 3) is used. Free-slip
conditions are applied at the lateral and upperboundaries, and convective outflow conditions areapplied at the outflow boundary. At the surfaces of theground and the isolated hill, no-slip conditions areimposed. The time step is set to ⊿t = 10-5 h/Uin (refer toTable 1).
In contrast to RIAM-COMPACT®, Meteodyn WTadopts the RANS equation with a one-equation closure
scheme (k-L turbulence model; here, k and L refer toturbulence energy and the turbulence length scale,respectively.) Details of the computational techniques
usedinMeteodyn
WT
areexplainedin
http://meteodyn.com.

zFlow

Inflow boundary

94

Characteristic velocity scale Uin (m/s)
Characteristic length scale h (m)(= elevation difference = z max - zmin)

Maximum elevation zmax (m)Minimum elevation zmin (m)Fig.4

Characteristic scales used in RIAM-COMPACT®

Reports of Research Institute for Applied Mechanics, Kyushu University No.149 September 2015

Flow

A flow separation region (vortexregion) characterized by negative
wind velocities does not form.

Fig.5

Streamwise (x) wind velocity distribution at the center of the span (y = 0),
Meteodyn WT, k-L RANS model, Re = 107

Flow

Fig.6

Streamwise (x) turbulence intensity distribution at the center of the span (y = 0),Meteodyn WT, k-L RANS model, Re = 107

A flow separation region (vortexregion) characterized by negative
wind velocities does not form.

FlowFig.7

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

95

96

Uchida et al. : Analysis of the Airflow Field with Commercially Available CFD Software
Flow(a) Instantaneous field(b) Time-averaged field
Note: evaluated for t =100 – 200 h/Uin

Fig.8

A flow separation region (vortexregion) characterized bynegative wind velocities forms.

A flow separation region (vortexregion) characterized bynegative wind velocities forms.

Streamwise (x) wind velocity distribution at the center of the span (y = 0) in the vicinity of the isolated hill,
RIAM-COMPACT®, the standard Smagorinsky LES model, Re = 106

3. Simulation results from the two software
packages

velocities clearly exists downstream of the isolated hillin the simulated flow at Re (=Uinh/ν) = 106.

Figs. 5 - 7 show results from the simulation with theMeteodyn WT software package (turbulence model: k-LRANS). These results (for a flow at Re = 107) indicatethat 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. Fig. 8 shows the resultsfrom the simulation with the RIAM-COMPACT®
natural terrain version software package (turbulencemodel: the standard Smagorinsky LES). Examinationsof the results reveal that a reverse flow region (vortexregion) characterized by negative values of wind

4. Summary
The present study compared the prediction accuracyof two CFD software packages for simulating a flow
around a three-dimensional, isolated hill with a steep
slope: 1) Meteodyn WT (turbulence model: k-L RANS),
which is one of the leading commercially available
CFD software packages in the wind power industry and2) RIAM-COMPACT® (turbulence model: the standard
Smagorinsky LES), which has been developed by thelead author of the present paper. Although the Reynolds
number settings differed between the simulation with

Reports of Research Institute for Applied Mechanics, Kyushu University No.149 September 2015
Meteodyn WT and that with RIAM-COMPACT®,
distinct differences in the flow patterns were identified
in the vicinity of the isolated hill (especiallydownstream of the isolated hill) between the flows
simulated by the two CFD software packages.

Acknowledgements
The present study was partly funded by the Joint
research with Hitachi Power Solutions Co., Ltd. I wish to
express my gratitude to the funding agency.

References
1)

2)

Takanori UCHIDA, Validation Testing of thePrediction Accuracy of the Numerical Wind Synopsis
Prediction Technique RIAM-COMPACT for the Caseof the Bolund Experiment –Comparison against aWind-Tunnel Experiment–, Reports of RIAM,Kyushu University, No.147, pp.7-14, 2014
http://meteodyn.com/Addendum
As discussed in the main text of the present paper, noreverse flow region (vortex region) characterized bynegative wind velocities was identified downstream of
the isolated hill in the result from the simulation withMeteodyn WT (k-L RANS model), in which theReynolds number was set to Re = 107 (the default value
for Meteodyn WT). Accordingly, in order to investigate
if a flow pattern similar to that simulated with MeteodynWT (k-L RANS model) can also be formed with the useof RIAM-COMPACT® (the standard Smagorinsky LESmodel), a simulation is conducted with RIAMCOMPACT® in which the method for setting thesurface boundary conditions is modified.
Specifically, the values of the three components of thewind velocity from k = 2 on the computational grid ofMeteodyn WT (k-L RANS model) (k = 1 corresponds tothe surface of the ground or the isolated hill in
Meteodyn WT) are assigned as the values for the surface
boundary conditions (k = 1) for the simulation withRIAM-COMPACT® (the standard Smagorinsky LESmodel) (Fig. 9). That is, this simulation is one in which
non-zero values are assigned for the wind velocity at thesurfaces of the ground and the isolated hill at allcomputational steps (Dirichlet boundary condition). The
Reynolds number, based on the height of the isolatedhill, is set to Re (= Uinh /ν) = 104, and the computational
grids shown in Figs. 1 and 2 are used. For the inflow
profile, a uniform flow profile is adopted in which the97

wind velocity does not change in the vertical direction.
The other boundary conditions are set using the same
methods as those discussed in the main text of thepresent paper. The time step is set to ⊿t = 2×10-3 h/Uin.
The results obtained from this simulation are comparedto those from a simulation which is identical to thesimulation described in this addendum, except thatno-slip conditions are applied at the surfaces of theground and the isolated hill, that is, the threecomponents of the wind velocity are all set to zero as
Dirichlet boundary conditions. For convenience, thesimulation in which non-zero wind velocities are applied
as a Dirichlet boundary condition is denoted as Case 1,and the simulation in which all three components of thewind velocity are set to zero is denoted as Case 2. A
comparison of the simulation results is shown in Fig. 10.
An examination of Fig. 10 reveals that no flow
separation region (vortex region) characterized bynegative wind velocities forms downstream of theisolated hill in Case 1, in which non-zero windvelocities are applied as a Dirichlet boundary condition.
In this case, a flow pattern resembling potential flowexists downstream of the isolated hill instead. In
contrast, in Case 2, in which all three components of thewind velocity are set to zero as a Dirichlet boundary
condition, a distinct flow separation region (vortexregion) forms downstream of the isolated hill. Thus, theflow pattern which forms in the vicinity of the isolated
hill varies significantly according to the velocity
boundary conditions applied for the surfaces of theground and the isolated hill.

Fig.9 Velocity vectors at the surfaces of the groundand the isolated hill (k = 1), simulation for which thewind velocity from k = 2 in Meteodyn WT is used as aDirichlet boundary condition, Case 1

98Uchida et al. : Analysis of the Airflow Field with Commercially Available CFD Software

A flow separation region (vortexregion) characterized by negative
wind velocities does not form.

A flow separation region (vortexregion) characterized by negative
wind velocities forms.

Uniform flow(a) Case 1: Simulation result from the casein which non-zero wind velocities are applied
as a Dirichlet boundary condition

Fig.10(b) Case 2: Simulation result from the casein which all three wind-velocity components are
set to zero as a Dirichlet boundary condition

Streamwise (x) wind velocity distribution in the vicinity of the isolated hill at the center of the span (y = 0),
RIAM-COMPACT®, the standard Smagorinsky LES model, Re = 104