Validation of the simpleFoam ( RANS ) solver for the atmospheric boundary layer in complex terrain

We validate the simpleFoam (RANS) solver in OpenFOAM (version 2.1.1) for simulating neutral atmospheric boundary layer flows in complex terrain. Initial and boundary conditions are given using Richards and Hoxey proposal [1]. In order to obtain stable simulation of the ABL, modified wall functions are used to set the near-wall boundary conditions, following Blocken et al remedial measures [2]. A structured grid is generated with the new library terrainBlockMesher [3, 4], based on OpenFOAM’s blockMesh native mesher. The new tool is capable of adding orographic features and the forest canopy. Additionally, the mesh can be refined in regions with complex orography. We study both the classical benchmark case of Askervein hill [5] and the more recent Bolund island data set [6]. Our purpose is two-folded: to validate the performance of OpenFOAM steady state solvers, and the suitability of the new meshing tool to generate high quality structured meshes, which will be used in the future for performing more computationally intensive LES simulations in complex terrain.


Introduction
Site selection for wind energy projects is a critical step of the decision making process for possible investment in wind energy.On-site measurement campaigns over areas of several kilometers square for wind farm siting can be very costly, and can only provide a local assessment of wind conditions.To overcome this hurdle it is currently very common to recur to computational fluid dynamics (CFD) simulations that can extend meteorological mast point source measurement data over an area of several kilometers.
Nowadays computational fluid dynamics (CFD) modeling for the atmospheric boundary layer (ABL) in complex terrain is becoming more ubiquitous, although it is quite time demanding to perform numerical solutions of the governing fluid equations without making simplifications to the flow dynamics.The most common approach used for ABL flows are the Reynold Averaged Navier ITM Web of Conferences Stokes (RANS) equations embedded with two-equation (k-) turbulence models [7], considering the flow incompressible, steady-state and turbulent assuming a neutrally stratified ABL.
In this contribution, we validate the CFD tool OpenFOAM (version 2.1.1)for ABL flows in complex terrain.We compare our simulations with the classical benchmark cases of Askervein [5] and Bolund [6] hills.Additionally, we test a new OpenFOAM-based tool for generating structured meshes in complex terrain [3].
The paper is divided in 4 parts.In Section 2 we describe briefly the two data sets studied.In Section 3 we discuss the equations and the numerical method used to model the wind field in the complex orographies.We discuss the results of the simulations in Section 4 and state some conclusions in Section 5.

Case studies
We study two classic benchmark tests for complex terrain: Askervein and Bolund hills.Askervein hill is a 116 m high (126 m above sea level), elliptical hill with a 2 km major axis and 1 km minor axis on the west coast of the island of South Uist, in the Outer Hebrides of Scotland.In 1982 and 1983 about 50 met masts for wind measurements were deployed for a measuring campaign [5].The measuring towers were placed in three linear arrays, one parallel and two perpendicular to the major axis of the hill.The field measurement points hill top (HT), central point (CP), together with the array deployment at 10 m height above the surface of the hill (lines A, AA and B) are shown in Figure 1 (top).The roughness length of the hill is z 0 = 0.03 m, and assumed constant all over the terrain in this investigation.
The Bolund experiment was a measuring campaign performed in 2007 and 2008.The Bolund hill, located north of the RISØ National Laboratory (Denmark), is a steep escarped island of about 130 m lenth, 12 m height and 75 m width [6].The island is connected to the coast by a narrow strip of land.The roughness length of the land is z 0 = 0.015 m, while the sea roughness length is z 0 = 0.0003 m.The two main measuring lines (line A and B), as well as the M0 reference mast position (for western wind directions) and M3 mast (top of the hill) are shown in Figure 1 (bottom).In this investigation, the narrow strip of land connecting the island to the coast is ignored.

Physics setup and solution method
The equations solved are the steady state RANS equations for the wind field [8].The topography of Askervein and Bolund hills is based on data provided by the Wakebench project [9].
The equations are discretized in computational domains of (streamwise × spanwise × height) 3 × 5 × 1 km 3 for Askervein hill and 300 × 200 × 300 m 3 for Bolund hill.A structured grid is generated with the new library terrainBlockMesher [3,4], based on OpenFOAM's blockMesh native mesher.We use about 5,400,000 hexahedral cells for the simulations presented in this paper.Additionally, we compare the new mesher with snappyHexMesh, OpenFOAM's unstructured mesher, with a maximum of 11,000,000 cells.
We use OpenFOAM's simpleFoam RANS solver.Two-equations turbulent models based on the kmodel are used for calculating the turbulent viscosity [7].
We use OpenFOAM's ABL libraries for setting up the boundary conditions.They follow the Richards and Hoxey (RH93, [1]) solution for the k-model at the inlet of the computational domain for the velocity U, turbulent kinetic energy k, and turbulent dissipation rate for a neutral ABL.This completely ignores thermal stratification.Coriolis effects are also ignored.
It is well known that the RH93 solutions are not entirely consistent with the k-model (see for example [2]).Recently, Gorlé et al [10] and Parente et al [11] proposed a set of boundary conditions which are fully compatible with the k-model solutions as proposed in RH93.These modifications include the addition of new source terms in the k-equations to make the equilibrium RANS equations match the RH93 solutions.In this paper, we choose to test different versions of the k-model available The boundary conditions are the following: • Incoming flow (Inlet boundary): Dirichlet conditions for velocity U, k, and using the Richardson and Hoxey expressions (ie, fully homogeneous fully developed 2D profiles), using the OpenFOAM libraries atmBoundaryLayerInletVelocity and atmBoundary-LayerInletEpsilon.For pressure, a Neumann zero gradient condition is assumed.• Outlet boundary.Flow is considered fully developed there, applying Neumann zero gradient condition for all variables, except for pressure.For this last variable a fixed (pressure = 0) value is set.• Sky (top boundary).Dirichlet conditions for U, k, and , prescribing the values corresponding to the inlet profiles at the height of the computational domain.Additionally, one can also use slip conditions here: zero for the normal component of a vector, and zero gradient for tangential, zero gradient for any scalar.• Ground (landscape terrain).Since we are only interested in modelling the inner region or constant stress layer, close to the viscous layer the variables are modelled via wall functions. 01002-p.3

ITM Web of Conferences
We use standard wall functions for k, and turbulent viscosity ν t [2].For U, no slip is assumed, with a Neumann zero gradient condition for p. • Sides (parallel to the flow direction).Slip conditions are used on these boundaries, which are assumed oriented parallel to the direction of the flow.
It is important to mention that we use standard wall functions for rough walls, with u = u * ln(Ez + /C s k + s )/ , where u * is the friction velocity, is the von Karman constant, E is an integration constant (E = 9.8), and z + is the non-dimensional distance from the wall, defined as z + = zu * /ν (with ν the kinematic viscosity).The wall function depends on the dimensionless roughness height k + s = k s u * /ν, where k s is the size of the sand grain.Following [2], the constant C s is calculated from a first order matching between the law of the wall and the inlet profile to ensure continuity for the wall functions and its first derivative, giving k s = Ez 0 /C s [12].

Askervein hill
We calculate the dimensionless speedup ratio, defined as (U − U inlet )/U inlet where, U is the horizontal velocity at a height of 10 meters above the local terrain along the line of masts and U inlet is a reference velocity at the inlet of the domain.Additionally, we calculate the normalized turbulent kinetic energy (TKE) as TKE/U 2 inlet .Figure 2a shows that different turbulence models have an influence in capturing the solution in the downwind side of the hill, in particular the k-and k-SST models.The RNG kperforms best in this region.As it can be observed from Figure 2b the simulation results using the k-SST model greatly underestimate the measured TKE data for both the upstream and downstream side of the hill.As for the k-model, simulation results indicate an overestimation of the measured TKE data in the upstream region, whereas the RNG k-model better predicts the measured TKE in this region.
In order to assess the quality of the structured mesh generated by terrainBlockMesher we do a comparison of the dimensionless speedup ratio along line A of Askervein hill, using both a unstructured snappyHexMesh mesh and a structured terrainBlockMesher mesh.Details of the mesh topology for both meshes are presented in Table 1.The results are presented in Figure 3a.All simulations use the RNG k-turbulence model.The snappyHexMesh mesh tends to overestimate the speedup at the downwind side of the hill, but the terrainBlockMesher mesh follows a closer trend with the measured data throughout the hill.All simulations underestimate the dimensionless speedup ratio at the top of the hill.Figure 3b shows the same comparison, but for the normalized TKE.The simulations using snappyHexMesh (4.5 and 11 million cells respectively) overestimate the TKE for upstream side of the hill.All simulations consistently underestimate the TKE 200 m from the hill (this last feature is commonly quite hard to capture, see [20] for more examples).However, the simulations using snappyHexMesh show a greater underestimation on the downwind side of the hill.

01002-p.4
First Symposium on OpenFOAM in Wind Energy   Finally, in order to check how sensitive are the results to the mesh size, we perform a mesh sensitivity study using the terrainBlockMesher with different horizontal mesh sizes.We only vary the horizontal resolution, keeping the same number of vertical cells.This is done in order to respect the recommendations for first cell location based on k s ∼ 20z 0 , as stated in [2], where k s is the sand-grain roughness height.We use the RNG k-model for the Askervein simulations.The results are presented ITM Web of Conferences in Figure 4, which show 3 different meshes: 0.5 × 10 6 , 1.5 × 10 6 , and 5.4 × 10 6 cells.For the speedup the curves barely differ from one another (Fig. 4a).There is a slight underestimation of the TKE in the downwind side of the hill for the 0.5 × 10 6 cells mesh, as can be seen in the closeup of this quantity shown in Figure 4b.

Bolund hill
In order to take into account the different roughness lenghts in Bolund island (z 0 = 0.015 m), surrounded by water with z 0 = 0.0003 m, we use an additional tool named identifyForestCells [3].This tool was originally developed for identifiyng cells which are to be modelled as a (porous) forest, using a mean flow model [21].In this case the tool is modified to select the shape of the island and generate a list of different roughness lengths which is then fed to the wall function for the turbulent viscosity ν t .
The same comparison presented for Askervein hill is presented here, using different turbulence models and meshing tools.We present the speedup factor, defined as (U − U 0 )/U 0 , where U is the horizontal velocity at a height of 2 meters above the local terrain along line B and U 0 is a reference velocity at the M0 mast position.Additionally, we show the normalized TKE, defined as (TKE − TKE ref )/U 2 0 , where TKE ref is a reference TKE at the M0 mast position.Note that we follow the same nondimensionalization as used in [22].
Figure 5 shows a comparison of the 3 turbulence models (standard k-, RNG k-and the k-SST).The three models perform in a very similar way, with a slightly better result obtained using the k-SST model for the speedup (Fig. 5a).On the other hand, the same turbulence model gives the worst results for the normalized TKE (Fig. 5b).All turbulence models underestimate the jump in TKE at the top of the hill.
Figure 6 shows a comparison between a snappyHexMesh and a terrainBlockMesher mesh with a similar number of cells (4.5 × 10 6 and 5.1 × 10 6 respectively), using the standard k-turbulence model.The snappyHexMesh simulation underestimates the speedup in most of the experimental points, and both mesh setups overestimate the data point at the upwind side of the hill (third experimental data point from the left in Fig. 6a).For the normalized TKE (Fig. 6b), the snappyHexMesh simulation follows more closely the experimental data, but both meshes underestimate the jump in the TKE in the 01002-p.6 First Symposium on OpenFOAM in Wind Energy  upwind side of the hill, with the underestimation being less pronounced for the terrainBlockMesher simulation.
A mesh sensitivity study is also performed for Bolund hill, using the terrainBlockMesher with different horizontal mesh sizes, much in the same way as it is done for Askervein in the previous section, keeping the same number of vertical cells and varying the horizontal resolution.We use the standard kmodel for the Bolund simulations.The results are presented in Fig. 7, which show 3 different meshes: 1.8 × 10 6 , 3.5 × 10 6 , and 5.1 × 10 6 cells.For the speedup the higher resolution of the 5.1 × 10 6 cells mesh follows more closely the experimental points (Figure 7a), the same being true for the TKE values (Fig. 7b).
Additionally, the residuals from the iterative solvers used in the solution of all the variables (Smoothsolver for the momentum and turbulence equations and Geometric-Algebraic Multi-Grid 01002-p.7 solver for pressure [23]) are presented in Figure 8.We show only the results for the coarsest grids: 0.5 × 10 6 cells for Askervein hill (Fig. 8a) and 1.8 × 10 6 cells for Bolund hill (Fig. 8b).All variables converged to residuals ∼ <10 −5 .

Conclusions
We have performed a validation of the simpleFoam RANS solver in OpenFOAM 2.1.1.Additionally, we have tested the new library for structured meshes terrainBlockMesher [3,4], based on OpenFOAM's blockMesh native mesher.The results from Askervein hill are very encouraging, with the numerical simulations following quite closely the speedup and normalized TKE experimental curves.Best results where obtained using the RNG k-model, in particular for the speedup on the downwind side of the hill.The simulations with First Symposium on OpenFOAM in Wind Energy structured meshes generated by the mesh generator terrainBlockMesher produced gave consistently better results than the simulations based on the unstructured meshes produced by the snappyHexMesh mesher.
For Bolund hill, the results are more disappointing.The influence of the three different turbulent models tested was minimal.The structured mesh generated using terrainBlockMesher produced results closer to the experimental data, although overestimating the speedup and underestimating the TKE at the upwind side of the hill.On the turbulence model side, recent studies [24] have shown better results using modified constants for the RNG k-model.The model constants for the other turbulence models were taken most of the time at their default values.A more careful study to select more appropriate constants is needed.For example, the inability of the k-model to provide good estimates in both the Askervein and Bolund hills can be related to badly chosen model and wall function parameters.On the mesh side, it is possible that a higher and more careful refinement is needed, especially on the highly escarped upwind side of the hill.
This study represents a first take at the Askervein and Bolund data sets, on the two main wind directions (210 • and 270 • ) of these data sets, for purposes of calibration and testing of the OpenFOAM solvers.A more detailed study is currently under way [25].The fully structured nature of the mesh produced by terrainBlockMesher makes it suitable for use in large eddy simulations.This will be tested in a future investigation.

Figure 1 .
Figure 1.Topography of Askervein hill (top), showing the lines where measuring masts were placed: major axes of the hill (A and AA), minor axis (line B), as well as the hill top (HT) and CP measurement point [5].Topography of Bolund hill (bottom), showing the measurement lines A and B, and the positions of the masts M0 and M3 (top of the hill) [6].

Figure 2 .
Figure 2. (a) Speedup (U − U inlet )/U inlet versus position along line A and (b) normalized turbulent kinetic energy (TKE) TKE /U 2 inlet versus position along line A. Three different turbulence models are shown: standard k-(blue circles), k-RNG (red triangles) and k-SST (green diamonds).The solid black line shows the experimental data.

Figure 3 .
Figure 3. (a) Speedup (U − U inlet )/U inlet versus position along line A and (b) normalized turbulent kinetic energy (TKE) TKE /U 2 inlet versus position along line A. Three different meshes are shown: a 4.5 × 10 6 cells mesh done with snappyHexMesh (sHM, blue circles), a 11 × 10 6 cells mesh done with snappyHexMesh (sHM, red triangles), and a 5.4 × 10 6 cells made with terrainBlockMesher (tBM, green diamonds).The solid black line shows the experimental data.

Figure 5 .
Figure 5. (a) Speedup (U − U 0 )/U 0 versus position along line B and (b) normalized turbulent kinetic energy (TKE − TKE ref )/U 2 0 versus position along line B, at 2 m height.Three different turbulence models are shown: standard k-(blue curve), k-RNG (red curve) and k-SST (green curve).The experimental data is shown as black squares.

Figure 6 .
Figure 6.(a) Speedup (U − U 0 )/U 0 versus position along line B and (b) normalized turbulent kinetic energy (TKE − TKE ref )/U 2 0 versus position along line B, at 2 m height.Two different meshes are shown: a 5.12 × 10 6 cells mesh done with terrainBlockMesher (tBM, blue curve) and a 4.45 × 10 6 cells mesh done with snappyHexMesh (sHM, red curve).In both cases, the standard k-model was used.The experimental data is shown as black squares.