Modelling complex terrain effects for wind farm layout optimization

The flow over four analytical hill geometries was calculated by CFD RANS simulations. For each hill, the results were converted into numerical models that transform arbitrary undisturbed inflow profiles by rescaling the effect of the obstacle. The predictions of such models are compared to full CFD results, first for atmospheric boundary layer flow, and then for a single turbine wake in the presence of an isolated hill. The implementation of the models into the wind farm modelling software flapFOAM is reported, advancing their inclusion into a fully modular wind farm layout optimization routine.


Introduction
Complex terrain issues are of high importance for on-shore wind farm planers. Often the local wind conditions are signicantly inuenced by the topography. A steep hill may for example cause speed-up on the top, as well as ow separation and recirculation at the downwind ank (see eg. [1]). While linearized solutions for low hills have been available for a long time [24] and are included in available wind farm layout codes, for many realistic sites with complex orography the determination of the optimal layout remains a challenge.
Wind farm calculation and layout optimization codes like PARK/UPMWAKE [5] and WAsP [6] (DTU, Denmark), FarmFlow [7] (ECN, Netherlands), WindFarmer [8] (DNV GL, United Kingdom) and FLaP [9] (University of Oldenburg, Germany) rely on fast calculation methods for the power produced by an ensemble of interacting wind turbines. They commonly invoke models that separate dierent aspects of the versatile and interacting physics involved, thus reducing the complexity of the problem. These may be models for isolated wakes, wake overlap, the wind energy converters and background wind elds. Models for mildly complex terrain, wake meandering, roughness and stratication eects also exist, some based on advanced calculation techniques. Here we present the implementation of novel Computational Fluid Dynamics (CFD) based models for complex terrain eects into the recently developed software flapFOAM [10], following a CFD based phenomenological approach. The general idea is to combine pre-calculated CFD results in a sensible way, avoiding the costs for solving systems of dierential equations during the run-time of the software.
Complex terrain eects both the background elds and the turbine wake evolution. In other words, hills contribute additional wakes to the wind eld, and they deect and deform those that are generated by upstream rotors. In full CFD studies, the interaction of wakes and topography can be solved simultaneously, and therefore consistently. For modelling codes that are based on the wake superposition approach, such as flapFOAM, one could consider the results of such simulations with all wind turbines switched o as background elds. However, the wake transformation due to the terrain still needs to be modelled explicitly. Additionally, for wind roses with many sectors, the pre-calculation of CFD background elds for each site of interest may be too costly, and explicit wake modelling appears preferable.
Many models of ow over low hills refer to the work of Jackson and Hunt [2], who provided an analytical solution for a two dimensional hump with small curvature in the atmospheric boundary layer. Their ndings are based on a linearization of the velocity perturbations due to the hill and near-wall equilibrium eddy viscosity. Taylor et al. [3] and Hunt et al. [4] continued and extended this work by heuristic adjustments and further theoretical analysis, respectively.
In the context of flapFOAM, hill models have to be generalized such that they accept arbitrary inow wind proles. This is required, since locally any combination of the background eld velocities and wake decits may approach the elevation. It is a fundamental requirement that all complex terrain models are combinable with all other model choices within flapFOAM, since the inter-comparison of models is one of the main purposes of the code. We therefore seek a description of the ow deformation due to the presence of the hill that is suciently general.
In this work, the required data on the ow behavior over hills is extracted from RANS simulations, rather than from analytical and eective models. Notice that this is not due to limits of the approach, and the implementation of existing hill models from the literature into flapFOAM remains for future work. Since RANS simulations are able to capture ow separation and recirculation eects, within limits that depend on the selected turbulence model, such eects can in principle be carried over to the modelling software. The corresponding hill models are therefore not limited to low elevation or high aspect ratios, as it is the case for linearized models.
The article is organized as follows. In Section 2, RANS simulations of the ow over isolated hills are presented. Section 3 then denes hill models and describes how they can be deduced from the previously obtained simulation results. The implementation of hill models in flapFOAM is briey noted in Sec. 4, before Section 5 gives results from a comparison of flapFOAM and full CFD wind elds, including a discussion. Finally, the article concludes in Sec. 6.

Flow simulations
The eect of isolated hills on neutral ABL ow is studied by solving RANS equations with OpenFOAM (version 2.1.1), initially for at terrain with uniform roughness length z 0 = 5 cm, and subsequently in the presence of a specic hill geometry under otherwise identical conditions. For both cases the mesh was obtained using the open source tool terrainBlockMesher [11], which creates structured meshes for complex terrain applications and has an option for adding hill geometries. In the present study all hills are based on the analytical hill shape from the EPA wind tunnel experiments RUSHIL [12] and RUSVAL [13]. This shape was generalized by considering the radius a γ a function of the azimuthal angle γ, dened by an ellipse as shown in Fig. 1, with minor and major axes directions e x and e y . This ellipse is conned to the ground, such that each three dimensional hill is determined by four parameters; one angle α dening the axes orientation, the hill maximal height H, and the minor and major radii L and L , respectively. These parameters for four hills that were studied in this work are listed in Table 1. The geometrical shape of the hill is then obtained as follows. For a xed angle γ, the ground distance from the centre l γ is translated into the continuous parameter ξ γ , 0 ≤ ξ γ ≤ a γ , by (compare [14]) Height proles along the e x axis for the hills listed in Table 1. and the height h γ at that position follows from The resulting height proles along the e x axis for the four studied hills are shown in Fig. 2. A hill whose ground ellipse is centered at x c = (x c , y c , 0) induces the coordinate transformation between global coordinates x = (x, y, z) and the coordinates in the hill reference frame x = (x , y , z ), cf. Fig. 1, where h(x, y) is the shape function of the geometry and α the angle between e x and e x . Notice that in the hill reference frame the main ow direction is given by Re x . Throughout this work, the z-direction unit vector e z is orthogonal to the ground and the main ow in all simulations is along the x-direction unit vector e x . For all simulated hills, listed in Table 1, e x and e y coincide with e x and the y-direction unit vector e y , respectively.
Mesh independence was checked in the symmetry plane parallel to the minor axis for hill no. 4 from Table 1    In all simulations the standard k − turbulence model was applied, with parameters C µ = 0.033, C 1 = 1.176, C 2 = 1.92, C 3 = −0.33, σ k = 1 and σ = 1.3. The inlet boundary conditions for k and throughout this work corresponded to the Richards-Hoxey solution [15], at the ground standard wall functions were used. All elds of all simulations converged with residuals below 10 −5 .
The velocity eld that describes undisturbed ow is denoted as U 0 (x ) in the hill reference frame, yielding U 0 (x) = R −1 U 0 (x (x)) with respect to global coordinates. The coordinate dependence may for example arise from the inow prole or upstream hills, wind turbine rotors or other obstacles. Likewise, U (x ) and U(x) denote the ow elds in the presence of the geometry, with respect to the hill and global reference frames, respectively. We furthermore dene the normalized velocity decit as for |U 0 (x )| > 0, and zero otherwise. Notice that the normalization is dened with respect to the local undisturbed ow. The simulation results for the e x -component of ∆u (x ) for the four studied hill geometries, denoted as ∆u , are shown in Fig. 3 at six dierent positons along the ow axis. The inow velocity was chosen as 10 m s −1 at 123 m height. For hill nos. 1 and 3 strong recirculation regions are observed behind the obstacle, as expected for small aspect ratios L/H.
For comparison, the ow over the hills from Table 1 was additionally simulated in the presence of a uniform non-rotating actuator disk and inow u ref = 10 m s −1 at hub height h hub = h ref , rotor diameter D = 114 m, power coecient c p = 0.412 and thrust coecient c t = 0.64. The actuator disk was discretized with 220 cells in all meshes. For these simulations the turbulence model was equipped with additional dissipation near the rotor to enhance the resulting wake eect [16]. 3. Flow models 3.1. Multi-wake ow in at terrain Let U at 0 (x) denote the time-averaged velocity vector eld that describes ow over at terrain. We then assume that the eect of adding N wind turbines to the domain can be approximated where S denotes the wake overlap function and ∆U at wake,i (x) = U at wake,i (x) − U at 0 (x) the wake decit associated with turbine i, where U at wake,i (x) represents the ow eld including the wake eect. Both S and ∆U at wake,i (x) are understood as explicit functions that are available either analytically or numerically, depending on the selected models, cf. [10].
Common choices for S are the linear addition of decits, S(d 1 , . . . , d N ) = n d n , or the addition of decit squares under the square root, S(d 1 , . . . , d N ) = n d 2 n . A popular example for the wake decit model is the Jensen model [17], for which the decit component in streamwise , and the other decit components vanish. Here R is the rotor radius located at x = 0, c t is the thrust coecient at the eective undisturbed inow wind speed U 0 at the rotor, and the wake radius expansion parameter is chosen as k = 0.07.

Flow over isolated hills
In our approach the eect of the hill geometry on the ow is twofold. First, the shape of the wake decit is geometrically deformed. With respect to hill coordinates, indicated by primes, this is modelled by a wake transformation function G , Notice that this function, just like the wake decits ∆U at wake,i , depends on all relevant ow parameters, like turbulence intensity, surface roughness and other factors. These parameters are treated as implicit constants of the models, and their dependencies are suppressed throughout this work.
As an example we consider the deformation G = ∆U at wake,i (x , y , z +(1−f )h (x , y )), referred to as the`shift wake transformation model', where f is the shift constant between zero and one, and h (x , y ) = h(x (x), y (y)). For α = 0 this implies ∆U hill wake,i (x) = ∆U at wake,i (x, y, z−f h(x, y)), whence the undisturbed wake decit is shifted upwards according to the product of the shift constant and the height of the hill. The eect of this model is demonstrated in Fig. 4, for the Jensen wake model, f = 0.5 and hill no. 1 from Table 1.
The second eect due to the presence of the hill geometry is a transformation of the local velocity vector. This can be expressed in terms of the normalized ow decit ∆u (x ), dened in Eq. (4), where in the presence of wind turbine wakes the undisturbed ow eld is We now assume that the normalized decit with respect to this background of deformed wakes can be captured by a dimensionless vector valued function T (ξ, η, ζ), depending on dimensionless coordinates, Again, this function has implicit dependence on the geometrical data of the hill and the underlying ow conditions, such as ambient turbulence intensity and the roughness length. While the latter can be understood as model constants when generalizing T from one hill geometry to another, the underlying geometrical scaling of the eect may be a strong modelling assumption.  We refer to the choice of T by the term`hill eect model', since it aects the wind eld as a whole, in contrast to the`wake transformation models' discussed above. The choice T (ξ, η, ζ) = 0 represents the`trivial' hill eect model, which simply imposes the coordinate transformation (3) on the undisturbed velocity eld U 0 (x ), which in the case α = 0 yields U(x) = U 0 (x, y, z − h(x, y)), compare the upper panel of Fig. 4. Notice that this model does not adjust the orientation of the velocity vectors, and neither speed-up eects nor recirculation areas are captured.
In order to obtain more realistic hill eect models, we follow an approach based on the interpolation of wind prole data or the pre-calculated CFD simulations from Section 2. For this, the full three dimensional CFD simulation results are rescaled to the given hill parameters. Assume that the normalized decit eld is known for a given hill geometry h cfd (x, y) with height H cfd , lengths scales L cfd and L cfd and orientation angle α cfd = 0, and is denoted by ∆u cfd (x ). Then the choice T (ξ, η, ζ) = ∆u cfd ξL cfd , ηL cfd , ζH cfd , generalizes the normalized decits to other hill geometries with α = 0. Notice that the coordinate transformations that are underlying ξ, η, ζ and ∆u cfd (x ) refer to dierent hill shape functions h(x, y) and h cfd (x, y) in that case. Likewise, the diering normalizations imply a rescaling of the simulated hill eect from the background eld that was used for the simulation to the background eld of the form (7). The hill eect models which are constructed from the numerical normalized wake decits of the hills 1 to 4 from Table 1 are referred to as`relDefHill1' to`relDefHill4', respectively. A demonstration for the Jensen wake model and hill no. 1 for the hill eect model relDefHill1' is shown in the lower panel of Fig. 4. It superimposes a rotation of the local vectors, speed-up eects at the hill top and a recirculation area onto the wake model ow.

Implementation in flapFOAM
The software flapFOAM is a wind farm modelling tool box that is currently in development at Fraunhofer IWES for the purpose of fast wind farm calculations and layout optimization.  It calculates the local wind eld based on the concept of wake superposition, as formulated in Eq. (5). In addition to standard wake models like the Jensen, Frandsen, Larsen and Ainslie models it also contains wake models that interpolate between sets of pre-calculated single-turbine wake results from steady CFD simulations [10]. A rst comparison to experimental data from a laboratory scale wind turbine model can be found in [18]. Further validation is in progress. For this work, the calculation of the total wind velocity at a given point was extended to the background (7), using (4), for selectable hill eect models T , cf. (8), and wake transformation models G, cf. (6). Due to the modular structure of the code, this is compatible with all model choices of flapFOAM, eg., with all implemented wake models. 5. Results and discussion

ABL ow over isolated hills
The ow of the ABL prole described in Section 2, xed by the condition of wind speed 10 m s −1 at 123 m height and roughness length z 0 = 5 cm, over hill no. 1 from Table 1 is shown in Fig. 5. This conrms that the implementation works correctly, since the simulation results exactly coincide with the hill eect model`relDefHill1'. Not shown are similar results for hill nos. 2, 3 and 4, which conrm this nding. In fact, due to the dimensionless display of the results, all lines in these graphs are exactly identical with the ones shown in Fig. 5, except for the one showing the CFD RANS simulations.
Naturally, the choice of the hill eect model matters for the proles over a given hill geometry. For hills with similar characteristics, eg. hills nos. 2 and 4 which have identical aspect ratio L/H = 4 and laminar ow behavior, the maximal error that is induced by selecting one of the hill models`relDefHill2' or`relDefHill4' is small. However, the models that correspond to hills nos. 1 and 3, with aspect ratios 2 and 1, respectively, give the wrong description for these hills, and vice-versa. The reason is the qualitatively dierent ow physics; the presence of a recirculation region and ow separation has a strong eect on the downstream proles behind the hill. Figure 6 shows wind velocity proles for the hills nos. 1 and 4 in the presence of a wind turbine of hub height 123 m at upstream distance of 1000 m, using the CFD-induced wake model of flapFOAM [10], with settings as described in Section 2. The displayed normalized decits are  taken with respect to the RANS solution for at terrain in the absence of the wind turbine. The wake transformation model was`shift' with shift factor f = 0.5 for both hills and all hill models. Also in the presence of a wake the hill models`relDefHill1' and`relDefHill4' give a good approximation of the full CFD results for the corresponding hills. There is a slight shift to smaller values of the decit component in ow direction at the hill top, for heights below 2H in both cases, and also for the upstream prole at x /L = −1 for heights 1 < z /H < 2.5. For the downhill prole at x /L = 1 we observe a similarly small shift in the positive ow direction for hill no. 1 (upper panel) at comparable heights, and in the negative ow direction below z /H = 1 for hill no. 4 (lower panel).

Wake ow over isolated hills
These discrepancies are due to the limits of the modelling approach of rescaling the eects from one ow situation, ie., pure ABL ow over the hill, to a very dierent one, ie., the ow of a wind turbine wake over the obstacle. However, the limits of the modelling possibilities have not been explored in order to overcome the mismatch in these regions. For this work, we simply compared the results for shift factors f = 0, 0.5, 1 and found best agreement at the hill top for f = 0.5. Additionally, a wake widening function may be introduced as a wake transformation model. However, in the end the models are required to reproduce ow elds behind hills for arbitrary wake situations, and a model netuning to a particular wake is not recommendable. A comparison of the model predictions for dierent wakes is beyond the scope of this paper and remains for future studies.
Both the results for ABL and wake ow over isolated hills show a signicant sensitivity on the selected hill eect model, ie., the hill geometry that is used to map the undisturbed ow eld onto the eld including hill eects. The reason is the very dierent behaviour of the ow depending of the aspect ratio L/H of the hill. This fact has to be represented by the models, especially when applying them to hill geometries with unknown outcome. One simple approach is to select from a data base of CFD based hill eect models, based on the aspect ratio of the geometry of interest, and the surface roughness. The dependency on details of the hill shape, and the possibility to combine basic shapes and their eects, remains to be studied. Finally, the transition from isolated hills to fully complex terrain remains to be addressed and validated with data from real wind farms in complex terrain.

Conclusion
We presented a CFD based phenomenological approach for the inclusion of ow modications due to the presence of hills in the context of the wind farm modelling code flapFOAM. The main assumption behind the presented models is that the local relative transformation of a given ow eld, as it is induced by the hill geometry, can be re-interpreted for dierent inow situations. As the main interest lies in modelling complex wind farms, the latter are characterized by overlapping wakes in an athmospheric boundary layer background. In this proof-of-concept paper the problem was reduced to a single wind turbine wake and an isolated hill, and the generalization to complete wind farms is the crucial next step along the presented line of work. Four dierent hill geometries were studied, and the resulting hill models were compared to full CFD RANS simulations. For both ABL and wake ows over the hills it was found that the models give good approximations of the total velocity decits, with the reservation that the results are sensitive on the ow physics of the underlying hill geometry. The presence of recirculation areas at the downhill ank and ow seperation has major impact on the proles at downstream positions, and the hill models fail to reproduce the simulated results if they do not include such eects. An automated selection of hill models with or without ow seperation, depending on hill parameters and inow conditions, remains for future work.
Similarly to the CFD based wake models of flapFOAM, the CFD based hill models leave the solution of momentum and continuity equations, as well as turbulence modelling, to the CFD solver. The modelling enters in the generalization of the results, and their superposition. All calculations then reduce to reading and interpolating pre-calculated data, which is fast enough for running layout optimization iterations for large wind farms on a single desktop computer. We hope that the presented approach may be a step towards the inclusion of complex terrain eects into such optimizations, without limited applicability to small altitude changes at the site.