Comparison of OpenFOAM and EllipSys3D for neutral atmospheric flow over complex terrain

. The ﬂow solvers OpenFOAM and EllipSys3D are compared in the case of neutral atmospheric ﬂow over terrain using the test cases of Askervein and Bolund hills. Both solvers are run using the steady-state Reynolds-averaged Navier–Stokes k – (cid:15) turbulence model. One of the main modeling differences between the two solvers is the wall-function approach. The OpenFOAM v.1.7.1 uses a Nikuradse’s sand roughness model, while EllipSys3D uses a model based on the atmospheric roughness length. It is found that Nikuradse’s model introduces an error dependent on the near-wall cell height. To mitigate this error the near-wall cells should be at least 10 times larger than the surface roughness. It is nonetheless possible to obtain very similar results between EllipSys3D and OpenFOAM v.1.7.1. The more recent OpenFOAM v.2.2.1, which includes the atmospheric roughness length wall-function approach, has also been tested and compared to the results of OpenFOAM v.1.7.1 and EllipSys3D. The numerical results obtained using the same


Introduction
Wind resource assessment is of major importance to assess the economic viability of a wind farm project. The traditional tools used by the wind industry rely on linearized flow models and do not perform accurately in complex terrain . Thus, there is a growing interest from the wind industry to use full nonlinear computational fluid dynamics (CFD) solutions.
Several CFD solvers are currently available to the industry. Among them, OpenFOAM (OpenFOAM, 2016) has recently received a lot of interest from both the research community and the wind industry. Its popularity seems to be firstly caused by its gratuity and its open-source license, which is in strong contrast to most industry-oriented CFD solvers, and secondly by its good performance 1 , e.g., at the Bolund blind comparison (Berg et al., 2011;. Open-FOAM therefore directly offers accurate results at no licensing costs. Two questions are, however, unanswered: how easily and how fast OpenFOAM can produce accurate results. These two questions are of course relative to other flow solvers, but they are relevant and important for the wind industry, as they can represent substantial costs in terms of both necessary hardware and manpower required, due to potentially extensive meshing and computational times. In order to address these questions, OpenFOAM is compared with EllipSys3D (Sørensen, 1995;Michelsen, 1992). EllipSys3D is an inhouse CFD flow solver designed from the ground up for wind energy applications (e.g., atmospheric boundary layer flows, flow over complex terrain and wind turbine rotor computations) at DTU Wind Energy (formerly Risø National Laboratory; Sørensen, 1995) and DTU-MEK (Michelsen, 1992).
The comparison presented here is focused on the mesh requirements of OpenFOAM relative to EllipSys3D and how this can affect the accuracy of the results. Furthermore, the two CFD codes are benchmarked in speed on the computer cluster at DTU Wind Energy, utilizing two Intel Xeon ® nodes (2 × 12 CPU cores) connected with InfiniBand interconnect.
In order to keep the comparison as close as possible, the two flow solvers are run, to the best possible extent, using the same models and parameters. The main model difference between the two flow solvers is the wall function applied for modeling the effect of the ground roughness on the flow. El-lipSys3D (Sørensen et al., 2007a) uses a wall function based on the atmospheric roughness length (Richards and Hoxey, 1993), while OpenFOAM v.1.7.1 uses the Nikuradse sand roughness model (Nikuradse, 1950). In OpenFOAM v.2.1.1 a wall function based on the atmospheric roughness length (Richards and Hoxey, 1993) is also implemented. It should be noted that difference in wall-function modeling has a significant impact on the mesh requirement concerning the height of the first cell above the ground level.

Basic equations
Both EllipSys3D and OpenFOAM are based on the finite-volume solution of the Reynolds-averaged Navier-Stokes (RANS) equations. If molecular viscosity is neglected due to the high Reynolds number and the eddy-viscosity hypothesis of Boussinesq (1877) is followed, then the RANS equations can be written as where the Einstein summation notation is used. u i = (u; v, w) denotes the mean velocity vector; x i = (x, y, z) are axes of the coordinate system, with z being the vertical direction; p is the dynamic pressure; f i represents body forces; S ij is the strain rate tensor; and ν T is the eddy viscosity which needs to be modeled. The transient term of Eq.
(1) has been retained even though the equations in this work are solved in the steady-state mode. The classical two-equation high Reynolds number kmodel (Launder and Spalding, 1974) is utilized in both flow solvers to calculate the eddy viscosity. Transport equations for the turbulent kinetic energy, k, and its dissipation rate , used by kturbulence model are solved and have the following form: The originally proposed model constants by Launder and Spalding (1974) were established for industrial flows, while slightly different values are appropriate for atmospheric flows (Bechmann and Sørensen, 2010 Tables 1 and 4). Identical constants have been used in both EllipSys3D-and OpenFOAM-based calculations.

Boundary conditions
In order to model the effect of surface roughness and avoid resolving the laminar sublayer, it is common practice to use the wall function for atmospheric flows. One of the major differences between EllipSys3D and OpenFOAM v.1.7.1 is related to how the wall function is implemented. In the more recent OpenFOAM v.2.1.1 the same wall function as the one in EllipSys3D is included in the official release.

EllipSys3D
In EllipSys3D the traditional high Reynolds number equilibrium assumptions are used to derive the wall function (see Sørensen, 1995;Sørensen et al., 2007b;and Hackman, 1982, for details). These inherently neglect the laminar sublayer, which usually is less than a millimeter thick for atmospheric flows. The logarithmic equilibrium profiles used for the mean wind speed, u, and turbulent kinetic energy are where z 0 is the roughness length, κ = 0.40 is the von Kármán constant and u * is the friction velocity. As seen from Eq. (4), the wall is placed on top of the roughness elements (u = 0 for z = 0) and is consequently displaced by the roughness length. This has been done to avoid a minimum height restriction of the first computational cell In order not to abandon the momentum equation in the near-wall cell by simply fixing the velocity according to log law (Eq. 4), EllipSys3D instead implements the wall function through the wall shear stress, τ 0 . Based on Eqs. (4) and (5) the shear stress is calculated from τ 0 = ρ u w 1 u w 2 using where z is the distance from the bottom cell face to the cell center. The equation is abandoned in the near-wall cells; instead, is specified according to while the k equation is reduced to a balance between production (P k = ν T |S| 2 ) and dissipation rate and a zero-gradient boundary condition is set for k. The boundary condition for the velocity is the standard no-slip condition ensuring the velocity to be zero at the wall.

OpenFOAM
OpenFOAM v.1.7.1 -Nikuradse's equivalent sand roughness length model The wall function available in OpenFOAM v.1.7.1 is based on Nikuradse's equivalent sand roughness length, k s (Nikuradse, 1950), and can be used to model flows where the laminar sublayer is important. The OpenFOAM implementation of an atmospheric equilibrium boundary-layer condition, according to the expression for surface roughness of Nikuradse, reduces to the following equation for the mean wind speed 2 (Martinez, 2011;Blocken et al., 2007), where C s = 0.5 is a roughness constant, E = 9.79 is a smooth wall constant and k s = E C s z 0 = 19.58 z 0 has been used. Comparing with Eq. (4) it is seen that the roughness for OpenFOAM v.1.7.1 is placed on top of the wall (u = 0 for z = z 0 ). While this is physically the "correct" approach, it does set some constraints on the heights of the near-wall cells. Firstly, as argued in Blocken et al. (2007), it is not physically meaningful to have grid cells within the roughness height; therefore, the height of the near-wall cell center should be at least equal to or larger than the Nikuradse's 2 Assuming C s k + s 1, where k + s = u * k s ν and ν is a molecular viscosity. roughness length ( z ≥ k s ). Blocken et al. (2007) do, however, state that it is mathematically/numerically possible to have z ≤ k s , and comparison between measurements and a sand-grain-based velocity function fit (Sumer, 2007) shows that their agreement is very good when z ≥ 0.2 k s . Secondly, OpenFOAM has a restriction on the cell aspect ratio, i.e., the ratio between the longest and the shortest cell dimension. Meshes with cell aspect ratio larger than 1000 have a tendency to introduce numerical errors and to converge slower (Martinez, 2011) and it is therefore difficult to make shallow cells that resolve the flow close to the ground. As a compromise between the cell constraints and the need to resolve the near-wall flow, z/z 0 ≈ 10.0 was found to generally give the least numerical error for flat-terrain simulations (Martinez, 2011) and is used in the following.
Using the built-in routine (nutRoughWallFunction), OpenFOAM v.1.7.1 implements the wall function by specifying the eddy viscosity, which, for fully rough flows, can be written as Similar to EllipSys3D, the equation is abandoned in the near-wall cells. Using the built-in function (epsilonWallFunction), is set to while a zero-gradient condition is set on k (using the kqRWallFunction). A no-slip condition is set on the velocity vector.
OpenFOAM v.2.1.1 - Richards and Hoxey (1993)-based model The model is identical to the model implemented in Ellip-Sys3D. Regarding the OpenFOAM implementation, the only difference compared to Nikuradse's roughness length model is in the way turbulent viscosity is determined. Using the built-in routine (nutkAtmRoughWallFunction), OpenFOAM v.2.1.1 implements the wall function by specifying the eddy viscosity as All other boundary conditions used in connection with this wall model are identical to the ones described in the previous paragraph.

Solution techniques
In the simulations carried out in this work, in both flow solvers, the RANS equations are discretized using  the QUICK scheme (Leonard, 1979) (QUICKV, the vectorial version for OpenFOAM). The k and equations are discretized using the first-order Upwind Discretization Scheme (UDS). The pressure is solved using the SIMPLE algorithm (Patankar and Spalding, 1972) and is accelerated using a multigrid approach (GAMG for OpenFOAM).

HypGrid
EllipSys3D uses the mesh generator HypGrid (Sørensen, 1998), a hyperbolic mesh generator developed at DTU Wind Energy. HypGrid creates a 3-D structured hexahedron volume mesh using a hyperbolic marching scheme, based on orthogonality and cell volume from a 3-D terrain grid surface definition. It can produce meshes with cells of low nonorthogonality and low skewness. There is no constraint on cell aspect ratio.

SnappyHexMesh
SnappyHexMesh is a hexahedral unstructured mesh generation tool included in the distribution of OpenFOAM. Snappy-HexMesh can use a 3-D terrain surface and iteratively build a mesh upon it. Some options allow construction of several cell layers of a controllable height above the terrain surface. This feature makes it possible to have a refined mesh in the region of high velocity gradients, close to the ground. Several regions can be selected to be refined to a desired level, and this method is creating a mesh of a cell aspect ratio close to 1. This is done in order to ensure that OpenFOAM can solve the numerical problems with the highest efficiency and accuracy using meshes generated with SnappyHexMesh. The meshing tool SnappyHexMesh can be run in parallel on a computer cluster or a PC.  Taylor and Teunissen (1987).

Simulation inputs
The inputs used to simulate Askervein Hill are fitted based on the upstream mast measurements.
The complete set of all simulation inputs, in accordance with Taylor and Teunissen (1987), is found in Table 1.

Description of the mesh HypGrid
The Askervein map is used as input to the in-house tool "surfgrid". The surfgrid tool generated a circular (32.5 km in diameter) surface mesh, with a refined region of 2 km × 2 km (4 blocks) in the center, at the position where Askervein Hill is located (see Fig. 1b). From this surface mesh the HypGrid tool is used to hyperbolically march the surface into a 6.5 km high 3-D mesh. The final 3-D volume mesh is composed of 20 blocks of 64 × 64 × 64 : 5.2 million cells (see Fig. 2b).
Two meshes are created with this method, one with a first cell center height of 0.4 m = 13.3 z 0 (ideal for Open-FOAM v.1.7.1 -the HG1 grid), and one with a first cell center height equal to a half of the roughness length z 0 = 0.03 m (ideal for EllipSys3D -the HG2 grid). Two additional meshes (with the first cell center height of z = 0.83 z 0 and z = 1.5 z 0 ) were created for OpenFOAM v.   Table 2. tion could, however, be obtained on those grids only if the mapFields tool is used to, for example, interpolate the final solution obtained on the grid with the first cell center height at z = 13.3 z 0 . As the proper speed convergence tests could not be conducted on the mentioned grids, the numerical results based on them are not included here. The coarse HG0 grid (Table 2) is created exclusively to accommodate a proper grid sequencing procedure in the Open-FOAM HypGrid-based computations.

SnappyHexMesh
The surfgrid-generated surface grid used for HypGrid meshing is converted to STL file format and utilized to generate a ground boundary surface suitable for the SnappyHexMesh OpenFOAM utility. The SnappyHexMesh-generated surface mesh close-up view is shown in Fig. 1a. A rectangular domain of 11.03 × 11.03 × 3.10 km is discretized using the blockMesh utility (see Table 2), creating a background mesh with a resolution of 95.9 × 95.9 × 38.75 m in the x, y, and z directions, respectively. Only the SHM4 mesh (Table 2) has a background mesh resolution of 77.5 m in the z direction. The cross section of the SnappyHexMeshcreated (SHM4) volume grid can be seen in Fig. 2a. Indicators, locating positions of refinement boxes and surface layers used to generate meshes SHM(0-4), are presented in Fig. 3. Changing refinement levels in refinement boxes 1-3 together with the number of inserted surface layers are the controlling parameters used in the SnappyHexMesh grid creation process ( Table 2). The refinementSurface and resolveFeatureAngle parameters are used to control the surface refinement level relative to the background grid. The coarse SHM0 (Table 2) grid is only intended for use in connection with grid sequencing procedures. As can be seen in Table 2, grids with optimal position of the first near-ground cell heights, suitable for both Open-FOAM v.1.7.1 and v.2.1.1, could be created in this way.

Simulation results
The main results of the Askervein test case are presented in Figs. 4 and 5. Both figures show the results of the simulations compared with the cup anemometer measurements along the line A -see Taylor and Teunissen (1987). The basic simu- lations are carried out on HypGrid-based grids HG(1-2) and SnappyHexMesh-based grids SHM(1-4). In Figs. 4a and 5a EllipSys3D calculation on structured mesh HG1 and Open-FOAM (both v.1.7.1 and v.2.1.1) computations on the same mesh are presented. The lines in red, black and blue therefore represent the results of OpenFOAM and EllipSys3D calculations using the identical mesh.
Overall, the results of all simulations are very similar to the measurements in the region upstream of the hilltop, both in terms of speedup and turbulent kinetic energy (TKE).
After the hill-top however, results of the simulations start to differ from each other significantly, in terms of speedup and to a lesser extend in terms of TKE. In comparison with the measurements, the speedup results are relatively close to the EllipSys3D simulation, using the mesh with a first cell height equal to the surface roughness (HG2 grid -dotted magenta line - One should also note the very good agreement between OpenFOAM v.2.1.1 and EllipSys3D (the dashed black and blue curves), based on an identical wall-function model (Richards and Hoxey, 1993) and calculated on identical computational meshes (HG1), especially in terms of speedup ( Fig. 4a) but also in terms of TKE (Fig. 5a).

Simulation time
Default values (i.e., from simpleFoam tutorials), regarding basic OpenFOAM solver inputs, are used in all OpenFOAMbased calculations. Several attempts have been made to change/tweak some of (many) multigrid pressure solver -GAMG parameters, but they all resulted in prolonged computational times and in some cases led to a periodic rather than monotonic decay residual behavior. For extensive details about all input parameters considered, the interested reader is referred to Martinez (2011).
The computational process in EllipSys3D has been done both with the standard five-level grid sequencing 3 and without it. In the OpenFOAM case the direct grid sequencing pro- Table 2. Askervein Hill: overview of different control parameters used to generate the SnappyHexMesh-created grids. The "blockMesh" column shows the number of grid points in the x, y, and z directions in the background mesh. The "Refinement surfaces" column shows the minimum level of surface refinement relative to the blockMesh-created background grid (first number in the parentheses) and the maximum surface refinement level (second number in the parentheses), which is used if a cell intersection angle (angle between two adjacent cells) is greater than resolveFeatureAngle (number shown in column three). The "Ref. box" 1-3 columns show the local grid refinement level relative to the background grid. The following three columns show grid size before the addLayersControls option is used, the number of added surface layers, and a total height of grid zone corresponding to added surface layers. The last two columns show the number of grid points in added surface layers (total number of grid points in the entire mesh shown in the parentheses) and the ratio between the first cell center height and the roughness length. The definition of all grids, including the HypGrid-based ones, together with the corresponding grid sizes is also included. Positions of refinement boxes, surface layers and background grid are indicated in Fig. 3. SHM0-4 are the SnappyHexMesh-based grids and HG0-2 are the HypGrid-based ones.
cedure does not exist. Here, two grids -the first consisting of 32 × 32 × 64 : 1.3 million cells -HG0 grid (Table 2) has been used in connection with OpenFOAM calculations on HG1 grid and the second one -SHM0 (Table 2) is used in connection with SnappyHexMesh-based OpenFOAM calculations. Upon reaching a suitable convergence level on HG0 and SH0 grids, OpenFOAM's mapFields tool has then been used to interpolate the results to the fine HG(1-2) and SHM(1-4) grids. The computational process on the fine grids is then continued until the convergence criterion is met. Also, the OpenFOAM computations are carried out without using the aforementioned procedure, i.e., solving on the fine grid using the standard initial values for all variables. The obtained computational times are presented in Table 3.
In the comparison of EllipSys3D to OpenFOAM runs (Table 3), especially the fastest ones obtained on grids of similar size and location of the first near-ground computational cell (the EllipSys3D HG1 grid and OpenFOAM SHM3 grid), it can be observed that EllipSys3D is approximately a factor of 1.9-2.5 times faster in obtaining the numerical solution of the same level of accuracy.

Simulation inputs
The inputs used to simulate Bolund Hill test case are based on quantities proposed by . The values used are presented in Table 4.   Bolund Hill front side caused severe difficulties with regard to obtaining a suitable OpenFOAM 3D meshes, without problems in cell face orientation. The grid validity check conducted using the checkMesh OpenFOAM tool showed that bad cells could always be located in the mentioned area. As the surfgrid tool does not have an option to visualize the grid during the creation process or provide a possibility to smooth the created grids, Pointwise's Gridgen mesh generation software has thus been used for the purpose of smoothing the surfgrid-generated mesh 4 . A close-up view of the generated surface grid is presented in Fig. 6b. The 4 km in diameter surface mesh is centered on and refined around the Bolund Island position. The HypGrid tool is used to hyperbolically march the grid in the third dimension (up to a height of 1 km). The final volume grid (Fig. 7b) is comprised of 24 blocks of 64 × 64 × 64 cells, approximately 6.3 million grid points. Basically, three meshes are created in this way, one with the first cell center height of z = 0.1875 m, i.e., k s (z 0(0.0003) ) ≈ 0.006 m < z < k s (z 0(0.015) ) = 0.294 (OpenFOAM v.1.7.1 -the HG1 grid), a second one with the first cell center height of z = 0.0125 m, i.e., z 0(0.0003) < z < z 0(0.015) (OpenFOAM v.2.1.1the HG2 grid), and a third one with the first cell center height of z = 0.0005 m, i.e., z 0(0.0003) < z < z 0(0.015) (EllipSys3D -the HG3 grid).
As can be seen, some compromises on the position of the first grid point in the surface normal direction (especially in the OpenFOAM cases) had to be made due to the change in the surface roughness length throughout the computational Wind Energ. Sci., 1, 55-70, 2016 www.wind-energ-sci.net/1/55/2016/   Table 5. domain. Even though no formal model-based restriction exists regarding the position of the first cell center height in the Richards and Hoxey (1993) implementation of the wall function in OpenFOAM v.2.1.1, the previously discussed aspect ratio issue dictated that z = 0.0125 m. With this z value the max aspect ratio is approximately 7000, and the only error/warning reported by the checkMesh tool was regarding the cell aspect ratio. Diminishing the z value further introduces several other mesh related errors reported by the checkMesh tool and any attempt to obtain the numerical solution resulted in almost immediate divergence.

SnappyHexMesh
A similar procedure to that in the Askervein Hill test case, regarding the re-use of the surfgrid-created ground surface mesh utilized in the HypGrid meshing tool, is conducted in the Bolund Hill case. Upon the surface grid conversion to STL file format, a rectangular domain of 2.3 × 2.3 × 1.0 km is discretized using the blockMesh utility (see Table 5) creating a background mesh with resolution of 30.7 × 30.7 × 16.7 m in x, y, z directions, respectively. The final surface grid created using the SnappyHexMesh tool is presented in Fig. 6a. Analogously to the Askervein Hill case, the cross section of the SnappyHexMesh-created (SHM3) grid is shown in Fig. 7a and indicators locating positions of refinement boxes and surface layers used to generate meshes SHM(0-3) are presented in Fig. 8. Similarly, both adjustment of the refinement levels in refinement boxes 1-3 and the number of inserted surface layers are the controlling parameters used for the SnappyHexMesh grid creation (Table 5). The refinementSurface and resolveFeatureAngle parameters are used to control the surface refinement level relative to the background grid. They seem to play a much more important role in the Bolund Hill case than previously, apparently due to the sudden and abrupt change in the surface topology characterizing the Bolund Hill case. It should be noted that increasing the surface refinement level to more www.wind-energ-sci.net/1/55/2016/ Wind Energ. Sci., 1, 55-70, 2016 than level 5 (Table 5) led directly to an inability of Snappy-HexMesh to create valid surface layers. The coarse SHM0 and HG0 (Table 5) grids are again only intended for use in connection with grid sequencing procedures.

Simulation results
The results of the Bolund test case are presented in Figs. 9 and 10. The computations are conducted for the incoming wind flow direction of 270 • , and results are compared with measurements along the line B (for details see . The obtained results are compared to the results submitted by other participants of the Bolund blind comparison (not shown here) and found to be in close agreement with the submitted numerical results attained utilizing two-equation turbulence models.
As a roughness change characterizes the Bolund Hill test case, some considerations regarding the position of the first grid point in the surface normal direction do exist. The wallfunction closure used in EllipSys3D flow solver does not restrict the position of the first grid point in the surface normal direction, so the position chosen here corresponds to 3 × roughness length of the water (z 0 = 0.0003 m) and 1/15 × roughness length of the land (z 0 = 0.015 m). In the OpenFOAM v.1.7.1 case, based on investigations of Martinez (2011), the largest roughness length for land had to be used as a basis for the grid creation process in order to avoid problems with the limitations of the Nikuradse's sand roughness length closure ( z/z 0−Land = 12.5). This, however, places the first grid point at the inlet and whole region upstream of the Bolund Island position at a very large relative distance from the terrain ( z/z 0−Water = 625). This issue could potentially negatively influence the results upstream of the Bolund Island position.
For this reason two different 2-D flat terrain computations involving the grids where the position of the first grid point in the surface normal direction was varied from z/z 0−Water = 12.5 to z/z 0−Water = 625 were conducted. The obtained results (not shown here), however, appeared to be in a close agreement with each other. Also, comparison of OpenFOAM (both v.1.7.1 and v.2.1.1) and EllipSys3D results at the position of mast 7 (M7) -positioned just upstream of Bolund Island -shown in Fig. 11 indicates that the mentioned issue does not seem to have any negative influence on most of the presented OpenFOAM (both v.1.7.1 and v.2.1.1) results. As seen from Fig. 11, only OpenFOAM (both v.1.7.1 and v.2.1.1) results based on the SHM1 grid deviate in terms of TKE from the rest of the computations. This issue will be discussed further in the next section.
From Figs. 9 and 10 a general good agreement can be observed between results of OpenFOAM v.2.1.1, v.1.7.1 and EllipSys3D together with their relatively good correspondence with the measurements. Largest differences between results can be observed in TKE plots at 2 m a.g.l. (above ground level). This subject will also be discussed further in Sect. 4.2. Furthermore, in particular, results using Open-Wind Energ. Sci., 1, 55-70, 2016 www.wind-energ-sci.net/1/55/2016/ FOAM v.2.1.1 (HG2 grid) and EllipSys3D, based on an identical wall-function model, have very good agreement. It is also seen that all calculations have similar difficulties in agreements with the measurements in the area immediately after the sharp Bolund Hill front edge, as previously reported in .

Simulation time
The solver inputs, apart from the Bolund Hill-specific ones presented in Table 4, regarding both OpenFOAM and Ellip-Sys3D, are kept identical to the inputs previously used in the Askervein test case. It should be noted that stable convergence could not be obtained on SnappyHexMesh-based grids using the QUICKV discretization scheme. For that reason, a formally second-order linearUpwindV 5 scheme was used in those cases. Computations on the HypGridbased HG1 grid in OpenFOAM were redone using the linearUpwindV scheme in order to be able to access the speed of convergence and quality of the obtained results in a proper manner. The OpenFOAM v.2.2.1 runs on the HG2 grid are conducted using the PCG pressure solver also (results in parentheses in Table 6) instead of GAMG. The grid sequencing procedure is analogous to the one presented in Sect. 3.1.4, only here a grid corresponding to grid level 2 in EllipSys3D (every second point in all three directions is removed -the HG0 grid) has been separately created and utilized for the mesh sequencing procedure in the OpenFOAM HG(1-2) runs. The coarse SHM0 grid is created for the same purpose according to the specifications in Table 5 with regard to SnappyHexMesh-based OpenFOAM simulations. The obtained computational times are shown in Table 6. Comparing EllipSys3D to fastest OpenFOAM runs on grids of similar size (Table 6), it can be observed that Ellip-Sys3D, in obtaining the numerical solution of the same level of accuracy, is approximately a factor of 4-6 times faster when the grid sequencing procedure is turned on and a factor of 1.3-1.7 times slower when it is turned off.

Discussion
Firstly, it should be pointed out that several simulations using the Nikuradse's equivalent sand roughness length wall model have been rerun in the more recent OpenFOAM v.2.1.1. Also, OpenFOAM v.2.1.1 has been compiled using the optimized Intel icc compiler for the Intel Xeon ® (enclosed in the DTU Wind Energy cluster facility) CPU platform, and selected simulations were rerun in this environment as well. In both cases no noticeable differences in the computational times, compared to the ones presented here, were observed.

SnappyHexMesh
OpenFOAM's own SnappyHexMesh meshing tool seems to have reasonable grid generation applicability and flexibility for the investigated neutral atmospheric boundary layer (ABL) flow over complex terrain, although a number of problems were encountered during the meshing process.
In the Askervein Hill case it was possible to create several grids with good general and boundary layer resolution capabilities using it. Dedicated grids for both Nikuradse's equivalent sand roughness length wall model and the atmospheric roughness length wall-function approach could be made directly. Generating the surface layer, crucial for appropriate boundary layer description in the ABL simulations, was a quite difficult task. Only usable results were obtained by splitting the add surface layer grid generation process form the rest of the SnappyHexMesh mesh generation procedure and disabling all mesh quality checks during this phase. Otherwise, very strange results with several regions of missing surface layer parts were obtained. The grids created during the present work basically reflect more of a limit in what SnappyHexMesh is capable of handling regarding the addition of surface layers in the grid, rather than a carefully considered user-based specification of sizes and extent of different surface layer parameters.
In the Bolund Hill case, some of the mentioned surface layer generation problems were even more emphasized. The abrupt change in surface structure at the Bolund Hill front side created severe problems in the gener-Wind Energ. Sci., 1, 55-70, 2016 www.wind-energ-sci.net/1/55/2016/ ation of the surface layers, so extending, for example, the refinementSurfaces level parameter to more than 5, in order to better approximate the Bolund ground surface, was not possible in the current case. Basically, there is very little freedom in specification of many surface-layer-related parameters. Generally, it seems very difficult as a user to be in control of the surface layer creation process using the Snap-pyHexMesh tool, making it quite difficult for it to be used consistently in relation to grid generation processes relevant for ABL flows.

HypGrid
In contrast, HypGrid, a meshing tool developed deliberately to create low-skewness hyperbolic 3-D meshes based on complex surface topologies, appears to be able to cope with both investigated geometries without significant problems. EllipSys3D could easily handle all grids created by HypGrid, but some adjustments, described in Sect. 3.2.2, were necessary in order to make suitable grids for OpenFOAM runs.

Askervein Hill case
The results obtained with OpenFOAM for the Askervein Hill case show very good general agreement with EllipSys3D and cup anemometer measurements in terms of speedup curves presented in Fig. 4b. In particular, results using OpenFOAM v.2.1.1 SHM(2-4) grids, OpenFOAM v.1.7.1 SHM(3-4) grids and the EllipSys3D HG2 grid (Fig. 4a) seem to have a very good correspondence, both prior to and after the hilltop. The OpenFOAM (v.1.7.1 and v.2.1.1) calculation based on the HypGrid-created mesh (HG1) together with EllipSys3D calculation on the same grid appear to deviate significantly from the abovementioned cases and measurements on the lee side of the hill. Results based on the Open-FOAM v.1.7.1 SHM1 grid seem to be placed in between the two above-described sets of results. An important thing to note here is a very good correspondence between results of OpenFOAM v.2.1.1 and EllipSys3D, which are run on identical computational grids (HG1) using the same wall-function modeling approach 6 . Regarding the TKE plots presented in Fig. 5, a good correspondence between all computations (and partially measurements) can be observed on the front side of the hill. Only the EllipSys3D HG2 grid-based calculation deviates from the rest of the computations in the immediate vicinity of the hill top. A similar behavior regarding the TKE in this particular region has been observed in the OpenFOAM v.2.1.1 calculations on previously mentioned grids with cell center heights of the order of the roughness length -z = 0.83 z 0 and z = 1.5 z 0 (not shown here).
The general deviations between numerical results on the lee side of the hill are much larger, but the differences between all computations and measurements here appear much more significant (numerical findings underpredict measurements by more than 50 %) and dominant than the differences between the computations. In the obtained steady state solutions occurrence of the flow separation has not been detected. As intermittent flow separation seemed to occur during the observational period, and as separation increases the general turbulence levels, this can be a possible reason why currently used RANS turbulence model cannot predict levels of turbulent kinetic energy on the lee side of the hill accurately (Undheim et al., 2006). The RANS models in general are reported to have a substantial problem in predicting the measured turbulent kinetic energy levels correctly in the mentioned zone (Sørensen, 1995;Kim and Patel, 2000;Eidsvik, 2005). However, Castro et al. (2003) -using high-order accurate schemes and unsteady RANS (URANS) formulation -seem to better capture the measured turbulence properties whereas recent large eddy simulation (LES) studies (Chow and Street, 2009) and hybrid RANS-LES studies (Bechmann and Sorensen, 2011) do show some general, significant and promising improvements in this regard.
Considering OpenFOAM results from both Figs. 4 and 5, the OpenFOAM computations, based on SnappyHexMeshcreated grids SHM(2-4), seem to have the best agreement with measurements and EllipSys3D.

Bolund Hill case
The first observation regarding the Bolund Hill results presented in Figs. 9 and 10 is very close agreement between results of OpenFOAM v.2.1.1 (HG2 grid, dashed cyan line) and EllipSys3D (HG3 grid, red line), which both use the Richards and Hoxey (1993)-based wall-function approach. This close agreement indicates that both flow solvers are indeed generally quite capable of producing reliable CFD results.
The differences in wall-modeling approach appear to play a dominant role in flow predictions in the recirculation zone on the lee side of the hill. Here, the OpenFOAM v.1.7.1based calculations (HG1 (HG1 LU) and SHM2) seem to collapse to a single curve, while OpenFOAM v.2.1.1-based calculations (HG2, SHM(1,3) grids and EllipSys3D HG3 grid) seem to collapse to another curve at both 2 and 5 m a.g.l.
Regarding the TKE plots in Fig. 10, it is seen that all simulations at the position of mast 6 (M6) seem to have difficulty in producing the correct level of the turbulent kinetic energy, especially for the measured height of 2 m a.g.l. Considerable variation between the results can be seen in this particular region, where almost all OpenFOAM v.2.1.1 and EllipSys3D results predict higher peak values close to the M6 position than the OpenFOAM v.1.7.1-based calculations. The SHM1 grid-based calculation (now in OpenFOAM v.2.1.1) deviates again from general tendencies and is seen to produce almost identical results to the OpenFOAM v.1.7.1 calculation on the same grid.
The behavior of computed results on SHM1 grid therefore appears more to reflect the ability of the grid to "snap" the correct Bolund ground surface (refinementSurface level 3 is used in SHM1 grid vs. refinementSurface level 5 in SHM(2,3) grids) rather than the wall-modeling approach used in the computations. This inability of SHM1 grid to properly represent the Bolund ground surface is also visible in results at mast position M7 (Fig. 11), especially on the TKE plot subfigure. It should be noted that increasing refinementSurface level to a value higher than 5 resulted in an inability of SnappyHexMesh to create any valid surface layers, due to grid distortion of closely spaced almost perpendicular cell layers, underlining a general difficulty/inability of full user control over surface grid layer creation in SnappyHexMesh OpenFOAM utility. Therefore, the OpenFOAM v.2.1.1 SHM3-based results represent the highest peak value close to M6 position which can be produced with the current SnappyHexMesh-based setup. The Open-FOAM v.2.1.1-based HG2 computation is seen to easily extend the SHM3 mesh predicted peek value in this particular region.
Comparison of OpenFOAM v.1.7.1-based (HG1 and HG1 LU) results which use identical settings and grid, with the only difference being in the discretization scheme used (QUICKV vs. linearUpwindV), shows practically no difference in the obtained results, indicating that utilization of the linearUpwindV scheme for SnappyHexMesh-based computations does not seem to impair the quality of the computed results.

Askervein Hill case
Computational times from Table 3 show that a converged solution on a SnappyHexMesh-based grid SHM3 is between 3.5 and 6.5 (grid sequencing on) and 5 and 7 (grid sequencing off) times faster to obtained than results on a HypGrid-based grid (HG1) with similar number of grid points using both the OpenFOAM v.1.7.1 and v.2.1.1 solvers. A closer look at the residual curves showed that a single iteration takes roughly the same amount of time in both SnappyHexMesh-and HypGrid-based cases, but the number of iterations needed to reach the convergence level of 2 × 10 −4 is, as indicated in Table 3, much higher in the HypGrid-based mesh case. This indicates that a structured meshing tool like HypGrid might not be the most optimal choice for grid generation purposes in OpenFOAM.
Focusing now on SnappyHexMesh-based results, it can be seen that a speed increase close to a factor of 2 is gained by using the grid sequencing procedure. The same is almost true in the HypGrid-based cases: there is slightly lower speed gain in the OpenFOAM v.1.7.1 case and slightly higher gain in the OpenFOAM v.2.1.1 case. In the SnappyHexMesh SHM3 case the speed differences between OpenFOAM v.2.1.1 and OpenFOAM v.1.7.1 computations are almost negligible, indicating that differences in the wall-modeling approach used in the two calculations do not seem to affect the convergence speed process significantly. Also, considering the increase in the number of grid points -following cases SHM (1,3,4) OpenFOAM v.1.7.1 and SHM(2,3,4) OpenFOAM v.2.1.1 -a proportional reflection in the computational times from Table 3 can be observed, although the proportionality factor seems to vary in a nonlinear manner.
Compared to the Bolund Hill test case, the grid sequencing procedure in the EllipSys3D computations does not appear to influence the computational time considerably. The reasons for this probably lie in the general flow complexity (or rather lack of the same) of the Askervein Hill case, as it apparently prevents the solver from fully utilizing the advantage of a solution on the coarser grid level, compared to using the standard start guess in the solution procedure on the fine grid. Placing the first near-wall cell very close to the ground (the HG2 grid) is seen to have a slight influence on the computational times also (relative to the HG1 grid).
In comparison of the EllipSys3D and OpenFOAM runs, best done on grids of similar size and position of the first cell center above ground level (i.e., the HG1 and SHM3 grids), it is seen that EllipSys3D is approximately 1.9 (grid sequencing on) and 2.5 (grid sequencing off) times faster in obtaining the converged solution.

Bolund Hill case
In the Bolund Hill case 7 , a significant difference in computational times between SnappyHexMesh-and HypGridbased OpenFOAM cases is again observed (Table 6). When comparing results on grids of similar sizes and near-ground surface resolution capabilities, conducted using the same linearUpwindV discretization scheme, a speed difference of a factor of 1.8 (grid sequencing on), 4.7 (grid sequencing off) for OpenFOAM 1.7.1 (cases HG1 and SHM1) and a factor of 3.4 (grid sequencing on) and 9.8 (grid sequencing off) for OpenFOAM 2.1.1 (cases HG2 -fastest, PCG-based solution and SHM3) can be observed. This underlines once again a general OpenFOAM issue with grids comprised of high-aspect-ratio cells. It is also seen from Table 6 that the multigrid GAMG pressure solver is not functioning optimally and is significantly outperformed by the more conventional PCG solver in HypGrid-based OpenFOAM v.2.1.1 cases.
Furthermore, Table 6 indicates that converged solution using QUICKV scheme is slower to obtain than the corresponding solution using linearUpwindV discretization scheme, especially when grid sequencing procedure is turned off (of the order of 20-30 %).
Grid sequencing procedures seem to have a smaller positive influence on the SnappyHexMesh-based calculations in the Bolund Hill case (speedup factor of < 2) compared to the Askervein Hill case, while the opposite is true for the HypGrid-based calculations (speedup factor of approximately 3-5).
The EllipSys3D code seems to benefit considerably from its automatic grid sequencing procedure, as it speeds up the computations by more than a factor of 10.
From comparison of the EllipSys3D run with best directly comparable OpenFOAM run (v.2.1.1 SHM3 smallest z, identical wall-modeling approach) it can be seen that Ellip-Sys3D is approximately a factor of 6 times faster in obtaining the converged solution (grid sequencing on), but despite the fact that OpenFOAM SHM3 computation includes almost 20 % more grid points, it is faster (approximately 30 %) to obtain the solution on it than using the structured EllipSys3D solver with grid sequencing turned off.

Conclusion
In this work, the unstructured OpenFOAM flow solver is compared to the structured EllipSys3D flow solver on two test cases calculating neutral atmospheric boundary layer flow over complex terrain.
Two meshing tools are considered: the structured hyperbolic 3-D mesh generator HypGrid and OpenFOAM's own hexahedral mesh generator SnappyHexMesh. OpenFOAM was found to be able to successfully perform calculations on the HypGrid-created meshes in both considered test cases. SnappyHexMesh was also able to produce reasonable grids in both the Askervein Hill and Bolund Hill test cases. A very important parameter for computational grids in ABL flowsheight of the first near-ground cell -proved to be very difficult to directly control using the SnappyHexMesh tool reflected in its (in)ability to create surface layers. This issue makes it quite difficult to use SnappyHexMesh consistently in relation to grid generation processes relevant for ABL flows, so a tool with more direct user control over this crucial part of the grid generation process can be recommended.
In terms of accuracy, both flow solvers are found to perform equally well on the two test cases, regarding both the mean flow velocity and turbulence quantities. In particular, OpenFOAM v.2.1.1 and EllipSys3D calculations, performed on identical (Askervein Hill case) and very similar (Bolund Hill case) computational grids, using the same approach to wall-function modeling of ABL flows (Richards and Hoxey, 1993), were found to have a great mutual correspondence, underlining the fact that both flow solvers are quite capable of producing reliable numerical results.
However, a large discrepancy in the speed performance is found. A very large difference in calculation times is ob-tained between HypGrid-and SnappyHexMesh-based Open-FOAM calculations, indicating that a structured meshing tool, typically creating grids with high-aspect-ratio cells in ABL flows, causes the OpenFOAM solver to perform inefficiently. Results of the present work show that this performance issue can be partially addressed by introducing a grid sequencing procedure in the OpenFOAM runs. Generally, the grid sequencing procedure had a very positive effect on almost all OpenFOAM computations and can be highly recommended.
In the comparison of EllipSys3D, which utilizes the grid sequencing procedure by default, and OpenFOAM SnappyHexMesh-based calculation times, the structured El-lipSys3D solver is found to perform approximately 2-6 times faster on grids with similar properties.
Using a combination of hexahedral and polyhedral cells the SnappyHexMesh tool can produce suitable grids in many relevant ABL flow cases and bring the computational times close to the level of the structured EllipSys3D code, but the inability to fully control ground surface approximation and associated surface layer creation, as in the Bolund Hill case, can potentially impair the quality of the obtained SnappyHexMesh-based results. On the other hand, use of the structured HypGrid solver, where terrain description and surface layer creation can be fully controlled, proved to be accompanied with a high speed performance penalty in the OpenFOAM runs. Thus, an unstructured mesh generation tool where both surface approximation and accompanied surface layer creation is better controlled might be an optimal meshing tool for ABL flow calculations using the OpenFOAM solver.