A numerical study on the flow upstream of a wind turbine on complex terrain

. The interaction of a wind turbine with the upstream ﬂow-ﬁeld in complex and ﬂat terrain is studied using Reynolds-averaged Navier-Stokes (RANS) simulations with a two equation turbulence closure. The complex site modelled is Perdig˜ao (Portugal), where a turbine is located on one of two parallel running ridges. Simulating various wind directions with and without rotor, the impact of the rotor on the ﬂow-ﬁeld upstream is determined. This is compared and related to simulations with sheared and uniform inﬂow. The induction zones forming for these two inﬂows agree to such an extent, that shear could be interpreted as linear perturbation to the uniform inﬂow solution. However, for complex terrain this is not the case, as it is highly dependant on ﬂow features caused by the topography and their interaction with the rotor. Separation in the lee of the ridge plays a crucial role, as it dictates the wind turbine wake trajectory which in turn governs the orientation of the induction zone.


Introduction
The flow coming towards a wind turbine rotor is continuously decelerated by the rotor's thrust force acting on it.The thrust is in turn a result of the aerodynamic forces acting over the rotor blades.The area over which this effect acts is also referred to as the induction zone.Medici et al. [1] found through RANS computations and wind tunnel measurements, that the influence of the turbine extends up to 6 rotor radii (R) upstream.Field measurements with wind lidars confirmed the presence of the decelerating region in front of the rotor, though the exact extent of it is disputed with Asimakopoulos et al. [2] putting it 7R upstream and Slinger et al. [3] at 3R.The IEC standards for power performance measurements [4] on the other hand assume the turbine to have a negligible effect beyond 4R upstream.Therefore this value also acts as guideline for industry.
Existing literature tried to quantify the induction zone solely in flat terrain, thus we address a situation where the topography considerably impacts the flow-field upstream of the rotor.Using RANS simulations in connection with an actuator disc (AD) representation of the rotor, we assess the difference in the induction zone in complex topography from flat terrain and uniform inflow.As a representative test case the site of the New European Wind Atlas (NEWA) at Perdigão (Portugal) is taken -here a turbine is located on top of one of two parallel ridges as shown in figure 1.

Computational method 2.1. Numerical setup
Changing the topography around a wind turbine is similar to altering the inflow profile to a rotor.Therefore simulations are preformed for uniform and sheared inflow profiles as well as for several wind directions at the Perdigão site.As the induction zone is governed by the rotor thrust, the comparability of the different inflow scenarios is driven by applying exactly the same thrust coefficient C T in each simulation.This condition is met with an actuator disc representation of the rotor with prescribed forces.The normal force acting over a sectional area ∆A of the disc is given by Here the free-stream velocity V {∞,∆A} acting over the area ∆A is determined by extracting the velocity at ∆A from a rotor-less simulation as shown in Figure 3. Ultimately subtracting the w/o rotor with rotor Figure 3. Method for determining the axial force F {N,∆A} for each section of the actuator disc by extracting velocities from a rotor-less simulation.
rotor-less flow-field from the one with rotor V R − V isolates the influence of the induction zone ∆ V .Forces are calculated for different thrust coefficients C T , namely 0.36, 0.64 and 0.89, to investigate the sensitivity of the induction zone to this parameter.The simulations are performed by solving the incompressible steady-state Reynolds-averaged Navier-Stokes (RANS) equations under neutral stratification.By omitting Coriolis and assuming neutral stratification the wind field becomes Reynolds number independent.The employed turbulence model, a modified k − [5] formulation, captures turbulence at both, terrain and rotor scales.The rough logarithmic velocity profile models the near ground flows [6].The friction velocity is prescribed at the inlet boundary with z 0 = 0.1 and u = 0.4.Additionally the wind direction is specified at the Perdigão site.The wind directions are selected to cover 60 • sectors, divided by 15 • steps centred about the directions orthogonal to the rotor ridge, equivalent to 53 • and 233 • .

Flow solver
The in-house finite volume code EllipSys3D solves the incompressible RANS equations over a discretised block-structured domain [7][8][9] with collocated variables.Solving convective terms using the QUICK scheme [10] and the SIMPLE method [11] for the pressure-linked terms of the Navier-Stokes equations.A modified Rhie-Chow algorithm [12,13] avoids decoupling velocity and pressure in the presence of discrete body forces originating from an actuator disc (AD).

Actuator disc
In the computational domain the prescribed rotor forces are distributed over a permeable disc using an actuator disc model [14].The disc itself consists of a polar grid with 33 points in the radial and 180 points in the annual direction.The intersections between the disc and computational grid determine in which cell the forces are applied.

Numerical domains
Figure 4. Impression of the Perdigão computational mesh along the transect of the ridge, passing through the turbine location.

Uniform & sheared inflow
A box domain with side lengths of 25 radii (R) is created for both inflow cases.This minimises domain blockage (π/25 2 = 0.5%).For uniform inflow the actuator disc is located at the domain centre and surrounded by a finely meshed box with 2.5R side lengths.Using a grid spacing of R/16 for the inner mesh surrounding the rotor, has yielded sufficiently accurate results in previous studies [15][16][17] for structured box domains in uniform and sheared inflow.The same finely meshed box surrounds the rotor in sheared inflow.However the rotor centre is shifted to 78 m above ground, equivalent to the hub height of the turbine at Perdigão.The wall cell is set to 0.06 m and grows from a height of z 0 [7] hyperbolically in the vertical direction until reaching the equispaced fine box mesh at 16.5 m above ground.For both grids the mesh spacing increases hyperbolically towards the outer domain edges starting from the finely meshed zones.The front, sides and top boundaries of the domain fulfil the Dirichlet and the rear the Neumann condition.The bottom face for uniform inflow is also of a Dirichlet type, whereas it becomes a no-slip boundary with a roughness length of z 0 = 0.1 in sheared inflow.

Perdigão
The turbine location acts as origin to an O-type mesh with a radius of 17 km.This meshing methodology allows to specify variable wind directions without changing the grid.
A finely meshed rectangular box surrounds the actuator disc with 8R side length (see figure 4).The inner mesh resolution matches that of the uniform and sheared domains.From the centrally located fine mesh the grid grows hyperbolically outwards, whilst following the terrain surface.The domain height is set to 10 km and discretised with 129 points.The in-and outflow boundaries are applied circumferentially to the O-mesh.The outflow region is defined as the arc resulting from a 90 • sector centred about the inflow direction.The velocities at the top and the inlet boundaries of the domain are prescribed, whereas at the outflow the velocity has zero-gradient.The roughness length is deduced from lidar canopy measurements taken at the site.The surface roughness as well as the topography is smoothed towards the edges to arrive at reference far-field quantities.
Error evolution in V with grid level relative to the finest mesh with spacing ∆x.The error is shown for the two main wind directions and calculated for an area upstream of the rotor.

Mesh sensitivity
The mesh spacings defining the grid surrounding the rotor at Perdigão are equivalent to those in the sheared inflow grid.However, the area it covers is larger incorporating a box with 8R side lengths.The area of interest studied in this paper extends 5R upstream and 1.2R radially outwards from the rotor centre.The discritisation error is calculated over this area for each grid level i with spacing 2 (i−1) ∆x, where ∆x represents the highest resolution and V the velocity magnitude: In complex terrain the error is highly depended on wind direction, due to the differing terrain interacting with the wind field.Therefore the error is shown for both main wind directions, 53 and 233 • , in figure 5. Convergence is achieved for both directions, though for 53 • the gain in accuracy from using the finest grid is larger.

The induction zone in complex inflow
In complex terrain the change in the flow-field caused by the rotor can be much more drastic depending on the effective topography.The latter changes markedly with wind direction as shown in figures 1 and 2.Even small differences in the topography can have a large effect on the induction zone, which is exemplary shown for wind directions 218 • and 263 • with C T = 0.89.For these directions the terrain gradient is larger leading up to the turbine than towards the other side.In figure 7 each row presents the results for one of these directions.The first two columns show the axial velocity component for undisturbed and rotor influenced flow-fields, respectively.In the final two columns contours of the normalised change in the axial velocity component are presented.Comparing the undisturbed flows, the much larger separation region downstream of the ridge for a wind direction of 218 • is noticeable.The separation is marked in light blue.Interestingly the transects for each direction are very similar, though.Once the rotor is introduced into the flow-field the separation bubble is pushed downwards, such that it almost disappears from the area shown for 263 • .Furthermore the acceleration at the hill top as well as its extremes are shown for the sectors centred around 53 • and 233 • .Furthermore the transect for each extreme is shown as reference.For the sector centred around 53 • the tendency is similar across all directions and the induction zone points slightly downwards.For the other sector the location of the maximum deficit covers a greater range over all wind directions.On average the induction zone points upwards, however.Apparently the induction zone does not follow the terrain and its direction is insensitive to changes in thrust.It is in fact the interaction of the rotor wake with the ridge flow that governs the shape of the induction zone.For a wind direction of 218 • a large separation bubble forms downstream of the ridge, which pushes the wake upwards which in turn leads to an upwards pointing induction zone.The correspondence of wake and induction zone directions agrees with the results of Brandlard et al. [18], who found the same behaviour for yawed rotors.Determining the induction zone in complex terrain is therefore immensely challenging, as it relies on predicting separation in complex topography correctly.The lower does not, but excludes some data lying far beyond the color range.Unsurprisingly the largest differences are mostly concentrated towards the lower part of the flow-field, close to the ground and towards the rotor edges.In case of sheared inflow the difference to uniform is remarkably low.This would most likely diminish further in-line with the roughness length.The equivalence between the uniform and sheared induction zone has been hinted at by Simley et al. [19].They measured the flow upstream of a full-scale wind turbine with a triple-lidar system in vertical and horizontal planes.The induction zones of both planes showed great similarity, when normalised by the hub-height wind speed or the shear profile, respectively.In complex inflow the differences are amplified compared to the sheared scenario.Nevertheless for a wind direction of 53 • the contours are clearly related to those for sheared inflow.This is not the case for 218 • .The large difference to uniform inflow arises from the upwards pointing induction zone, that ultimately stems from the large separation region in the lee of the ridge.Figure 10 shows the dependency of the normalised absolute difference on wind direction.For each complex terrain simulation the mean and extrema of an area bounded by −2 < x/R < −0.25 and −0.85 < z < 0.85 are determined for C T = {0.36,0.89}.Note that the two wind directions depicted in figure 9 are representative of the extrema in the mean value.Interestingly there is no clear dependency on thrust and despite the shallower terrain upstream of the rotor for the wind sector centred around 53 • , their values are on par with those of the other sector.
Ultimately it seems that the induction zone acts as a linear perturbation to the flow-field in sheared inflow, whereas in complex inflow the interaction is highly non-linear.

Conclusion
To the author's knowledge this is the first study on the induction zone in complex terrain.By comparing the induction zone between uniform, flat and complex inflow the influence of topography on the induction zone is assessed.There is generally strong agreement between the wind turbine induction zones forming with uniform and sheared inflow, whilst in complex terrain the topographically caused flow features play a crucial role.The speed-up effect in conjunction with separation forming downstream of the ridge can influence the shape of the induction zone.The latter influences the wake trajectory, which ultimately determines the orientation of the induction zone.Neither Coriolis nor stratification different from neutral have been considered in these simulations.Both can have an extensive influence on the wind turbine wake and separation and thus on the induction zone.Nevertheless the results of this study serve as a first indication of the dominating factors governing the rotor induction zone in complex terrain.A future study

Figure 1 .
Figure 1.Height contours of the parallel ridges at Perdigão (Portugal).The wind turbine location (o) and the investigated wind direction sectors are marked.

Figure 2 .
Figure 2. Height contours of the immidiate area sourrounding the wind turbine (o).

Figure 6 .
Figure 6.Contours of the normalised difference in the axial velocity component between the undisturbed and rotor simulation for uniform and sheared inflow with C T = 0.89.The ground is marked with a balck line in the sheared contours.

Figure 7 .Figure 8 .
Figure 7. Contours in each row correspond to solutions for the axial velocity component for wind directions 218 • and 263 • , respectively, and C T = 0.89.From left to right: Undisturbed flow; flow-field with rotor; difference between the former two with log-scale, zoomed difference with black line indicating maximum deficit as function of x.

3. 3 .Figure 9 .
Figure 9. Contours of the normalised absolute difference in the induction zone relative to uniform inflow.

Figure 10 .
Figure 10.Dependency of the normalised absolute difference on wind direction.For each complex terrain simulation the mean and extrema of an area bounded by −2 < x/R < −0.25 and −0.85 < z < 0.85 are shown for differing values of C T .
The Science of Making Torque from Wind (TORQUE 2016) IOP Publishing Journal of Physics: Conference Series 753 (2016) 032041 doi:10.1088/1742-6596/753/3/032041should systematically investigate different commonly encountered topographies.Consequently the uncertainty in the induction zone model should be determined as a function of the terrain.