Wind Farm Simulation and Layout Optimization in Complex Terrain

This work reports on incorporating complex terrain into wind farm simulations for the purpose of layout optimization. Adding complex terrain boundary conditions to NREL’s medium fidelity Reynolds averaged Navier-Stokes flow model, WindSE, produces significant separation, flow curvature, and speedup effects that would otherwise be difficult to capture with lower-fidelity models or a flat-terrain assumption. These flow features, in turn, can significantly impact the optimal turbine array layout. We demonstrate the addition of complex terrain improves agreement with SCADA data by 34.8% in mean average error. Additionally, we discuss modifications to the code that enable gradient-based layout optimization using terrainaware adjoints with arbitrary topography. Through several optimization case studies, we show that the layout optimization process takes advantage of speedup effects on terrain high points, and leverages flow curvature effects that modify wake trajectories. This yields substantial power improvements over gridded layouts, and hints at future research directions in simulation and optimization for wake trajectories in complex terrain.


Introduction
When planning a new wind farm, one of the major consideration is how to optimize the turbine layout to maximize estimated power output. Traditionally, this optimization is performed using engineering flow models that simplify key turbulent flow physics to obtain reduced computational costs. The simplifications in these models, such as linearization and wake superposition, are often appropriate in flat terrain, but in complex terrain many of their assumptions regarding thin shear layer approximations and attached parallel flows are violated. Furthermore, many terrain-induced effects cause local accelerations in the flow that can be beneficially employed by a sufficiently high fidelity optimization. Consequently, developing optimization-oriented models with higher fidelity than traditional tools can lead to reduced uncertainty in power outputs and improved layout design in complex terrain.
To remedy these modeling limitations, many in the wind energy community have begun employing computational fluid dynamics (CFD) models such as Reynolds-averaged Navier Stokes (RANS) models or large eddy simulations (LES). These CFD models preserve more of the true flow physics by solving either the ensemble averaged (in the case of RANS) or the spatially filtered (in the case of LES) Navier Stokes equations. In this study we focus our attention on RANS models as they strike a balance between computational cost and fidelity while improving the prediction of time averaged flow fields and power production. Furthermore, RANS models are increasingly being used in wind farm optimization studies, particularly when gradients are available through the use of adjoint solvers. In this paper, we demonstrate a layout optimization in complex terrain using an open source CFD tool from the National Renewable Energy Laboratory (NREL) called WindSE. As an adjoint RANS code, WindSE was previously used for gradient-based layout optimization on 3D domains with the assumption of flat terrain [1], as well as data-driven turbulence modeling [2,3] and active subspace dimension reduction [4]. However, the modifications reported in this paper are able to capture terrain-induced effects such as flow separation and wake curvature. These effects have a noticeable impact on the performance of a given wind farm layout, and ultimately lead to optimal turbine placement that is different than under the assumption of flat terrain. In Section 2 we provide a background on the RANS modeling in WindSE, its mixing length closure, its adjoint capabilities for obtaining high dimensional gradients, and review modifications to RANS models for complex terrain that we employ and those reported elsewhere in the literature. In Section 3 we report on simulation results for a wind farm in complex terrain that demonstrate the terrain-induced wake effects. In Section 4 we construct a layout optimization case study that demonstrates the impact of these terrain effects on the layout optimization results in terms of turbine locations and power performance. Finally, in Section 5 we conclude with a discussion of the simulation and optimization results.

Background
RANS modeling has become an increasingly popular choice for wind farm simulations because it strikes an attractive balance between computational complexity and physical fidelity. In recent years, RANS modeling has received particular attention on its handling of complex terrain and high dimensional optimization. Adjoint optimization techniques have enabled the use of sample-efficient gradient-based optimization techniques even when running expensive CFD models. These have been applied to wind farm applications using automatic differentiation for the adjoint system [1], a continuous adjoint formulation [5], and for model predictive control [6]. For complex terrain, a number of recent experimental and computational studies have examined wind turbine wakes in complex terrain. These studies have demonstrated the use of many different RANS turbulence models such as k − and k − ω [7], shown that wakes from turbines on hills can recover faster due to enhanced momentum entrainment [8], that RANS actuator disk models can compare favorably to LES with actuator lines [9], and found generally good agreement between field experiments, wind tunnel experiments, and CFD [10]. Initial studies on layout optimization in complex terrain have involved a hybrid CFD/engineering flow model approach with a random search [11], annual energy production with a random search [12], and CFD models with mixed integer programming [13]. To our knowledge, this is the first layout optimization study performed in complex terrain using adjoints and gradient-based optimization techniques.

WindSE
WindSE uses a Reynolds-averaged Navier Stokes (RANS) turbulence model to simulate the wind flow within a given domain and turbine layout. This system of equations is constructed in Python using FEniCS as its finite element backend [14]. Additionally, dolfin-adjoint, another Python package, is leveraged to automatically compute derivatives with respect to control parameters [15,16]. This allows WindSE to accurately optimize power output for high dimensional control inputs such as turbine location, yaw, and axial induction. A recent addition to the WindSE package is the ability to represent complex terrain as bottom boundary conditions in the RANS and adjoint solvers. To further account for varying inflow speeds and directions as part of a complete annual energy production (AEP) optimization, we additionally introduce modifications to the computational mesh and inflow boundary conditions to allow for variable inflow directions.

PDE Constrained Optimization
PDE Constrained Optimization consist of three major components: the objective function, the forward problem, and the adjoint problem. The objective function defines the quantity of interest that is being optimized. The forward problem consists of all the steps required to compute the current value of the objective function. The adjoint problem is used to find the gradient of the objective function with respect to the control variables. For WindSE, the objective function is the power output of the wind farm, the forward problem is the series of RANS simulations required to calculate the power, the control variables are the turbine coordinates, and the adjoint problem is derived by dolfin-adjoint.
WindSE solves the steady-state RANS equations with a mixing length closure for the Reynolds stresses and an actuator disk body force representing the wind turbines. For the studies reported here, we have assumed the turbines are operating in Region 2. Layout optimization is accomplished by maximizing the power output summed over a specified set of inflow directions and speeds, subject to the RANS equations as a PDE constraint as follows: whereê n is the unit normal to the rotor plane for the n th turbine, f AD,n is the turbine body force, A n is the turbine rotor swept area, c t,n and c p,n are the modified thrust and power coefficients defined for the local rotor disk velocity [17,18], ρ is the fluid density, k is the index over a set of inflow direction and speed conditions with corresponding weights α k that determine the probability of each inflow condition, and ϕ n (m n ) is a geometric smoothing kernel describing the actuator rotor disk defined as: where m n = [x n , y n , z n ] is the location of the nth turbine's hub, W is the width of the disk of influence, and R is the rotor radius. This method can be easily modified to account for arbitrary yaw.
The mixing length, mix , varies with height above ground according to the following relationship [19]: where κ = 0.41 is the von Karman constant, z is the vertical distance to the ground, and λ = 15m is the value of the mixing length in the free atmosphere [20,21].

Boundary Conditions
The RANS equations are solved using a log law inflow profile with a shear of 0.15 as the inflow boundary condition to the inflow plane and spanwise planes. A no-slip boundary condition is applied on the bottom surface, and a no-penetration boundary condition is applied on the top surface. A no-stress outflow boundary condition is applied at the outlet. To simulate different inflow directions we reuse the same computational mesh, but reassign boundary conditions to reflect changing inflow or outflow planes and velocity vectors.

Numerics
WindSE is built on top of the FEniCS project [14] and takes advantage of the dolfin-adjoint package [15,16]. The FEniCS project is a python based finite element toolkit that performs Just-In-Time compiling to C++, which combine to create an easy to use API with the speed of a C++ based code. The main advantages to using FEniCS as a backend, comes from its ability to assemble arbitrary PDE systems and its access to several solver methods. Currently, WindSE uses the MUltifrontal Massively Parallel sparse direct Solver (MUMPS) [22,23] to handle solving the forward problem combined with the Portable, Extensible Toolkit for Scientific Computation (PETSc) Scalable Nonlinear Equations Solvers (SNES) [24,25,26] to automatically handle the nonlinearity.
The optimization aspect of WindSE is handled by the dolfin-adjoint package, which was specifically built to directly interface with FEniCS. This package records a tape of every action required to perform a single forward solve, and then uses adjoint methods to estimate derivatives with respect to control variables such as turbine position, yaw angle, and axial induction factor. This information is then pass to another python package, scipy, which performs the optimization. The specific optimization is performed using Sequential Least SQuares Programming (SLSQP). This method allow for adding a minimum distance constraint that prevents two turbines from getting too close. For the studies in this paper, a constraint of two rotor diameters was used. Additionally, the location of the turbines is bounded by a farm constraint, which defines the maximum extents of the wind farm.

Complex Terrain Validation Study
We first investigate the simulation of an as-built wind farm in complex terrain without layout optimization. This simulation consists of using real terrain and turbine placements along with a log inflow profile with a velocity of 8 m/s at hub height blowing from the north. The turbines are model with a hub height of 80 m and rotor diameter of 80 m. Figure 1 shows the plant in consideration, with the relative terrain heights represented in the leftmost column. The column second from the left shows the streamwise velocities and the two rightmost columns show the spanwise and vertical velocities, with coloring to accentuate terrain-induced deviations from strictly streamwise flow. The bottom row show the results of the same simulation without terrain, highlighting how terrain affects the flow profiles. All velocities shown are measured from a constant distance of 80 m above the ground.
These results demonstrate that the terrain induces vertical and spanwise velocities that are not present in the flat simulation, resulting in changes to normal rotor inflow angles and wake deflections on the order of a few tens of meters. The terrain effects distort the normal wake shapes, transport the wakes vertically and laterally, and create local speedup effects on ridgetops. These effects are not found in engineering flow models but are critical to effective layout optimization and wake steering strategies.
In order to validate this simulation, Supervisory Control And Data Acquisition (SCADA) data was used for this site. First the SCADA data was filtered for conditions that matched the inflow conditions of this simulation. To do this, Turbine 7 from Figure 2 was chosen as a reference. Then the data was filtered for the time steps where the information for Turbine 7 closely match the simulated conditions with tolerances of ±0.25 m/s for the wind speed and ±15 • for the wind direction. Once these time steps were found, the data for the other 15 turbines was compiled. Finally, the filtered data was plotted in Figure 2 alongside the simulation data with and without terrain.
In Figure 2 the flat terrain simulation shows more or less constant power output. The only turbines that deviate are those that are directly waked by the turbines upwind. However, the simulation with terrain shows greatly improved agreement with the SCADA data. The mean average error (MAE) for the normalized power relative to the median SCADA data observations is 0.22 for the flat terrain simulations and 0.14 for the terrain-aware simulations, a reduction in MAE of 34.8%. The largest exceptions are Turbines 13 and 16 where the terrain model reduces error relative to the flat simulation but it still falls outside the SCADA data interquartile range. We attribute this discrepancy to two factors. First, there are turbines outside the computational domain whose wakes are reflected in the SCADA data but not captured in the simulation. Additionally, the simulation is configured to only allow horizontally homogeneous inflow conditions, but the inflow angles were observed to vary spatially in the SCADA data. Despite these limitations, the terrain-aware implementation of WindSE is able to significantly improve the agreement with field observations.

Layout Optimization Case Study
With a complex terrain simulation capability in place, we proceed to two optimization case studies. In the first case, we examine optimization of two turbines on a simple Gaussian hill to confirm that WindSE accurately computes the gradient of the power output with respect to  position. This requires WindSE and dolfin-adjoint to correctly handle the implicit change in elevation and corresponding increase or decrease in hub height windspeed caused by a lateral motion of the wind turbines. The second case study considers a layout optimization with four inflow directions and a skewed hill that breaks the axisymmetry of the first case study.

Case Study # 1
The first case study implements a simple Gaussian hill. This topography is chosen because it has a simple analytical expression, allowing us to compare results against an interpolated terrain representation. The hill is described by where (x 0 , y 0 ) is the center of the hill, and A is the maximum height of the hill. The coefficients a, b, c control the size and shape of the hill and are defined by: where, σ x and σ y control the extents of the hill in the x and y directions, and θ controls the skew. For this case study, the hill is constructed using the parameters in Table 1. WindSE has the ability to use the terrain data in analytical form shown in Equation 9 or to use an interpolation function to read in terrain data directly from file. The interpolation method allows for arbitrary terrain, and led to the same optimization results.  To verify the optimization capabilities of WindSE, two turbines are placed to the west of the hill located at (-520, 250) and (-520, -250). They have a hub height of 80 m and a rotor diameter of 80 m. Additionally, we optimize over two wind inflow directions, one from the west and another from the south with a wind speed of 8 m/s. Both inflows are assumed to be equally likely. With a single inflow direction, the optimal placement of the turbines would would simply be on top of the hill in a north-south line to take full advantage of the higher wind speeds. Adding the additional wind inflow that is offset by 90 • forces the optimal placement of the two turbines to be staggered to avoid one turbine being in the wake of another. After 18 iterations of WindSE's layout optimization this is exactly what we see shown in Figure 4.

Case Study #2
For the next case study, we consider a skewed Gaussian hill as defined by the parameters in Table 2. A contour plot of this hill is found in Figure 6. Additionally, this case study involves optimizing a layout of 9 turbines across four wind inflow directions. The turbines have a hub height of 80 m and a rotor diameter of 80 m. The inflow direction are simply from the north, south, east, and west. The inflow speed is 8 m/s. All inflow conditions are assumed to be equally likely. Future studies will expand on this simulation by including more turbines and adding additional controls such as yaw and axial induction. Additionally, the optimization is performed with the same setup but without terrain to highlight the optimization impacts.

Parameter Value Units
A 100 m x 0 0 m y 0 0 m σ x 200 m σ y 400 m θ π/4 rad Table 2. These are the parameters used to generate the skewed Gaussian hill.
The initial gridded layout and optimization results are shown in Figure 6. When the optimization is perform without terrain, there is a 49.0% improvement in power production over the gridded layout. When the skewed hill is considered, WindSE improved the objective function by 123.2% over the gridded layout. This was achieved by placing 5 of the 9 turbines on the top of the skewed hill and by ensuring minimal wake interactions by staggering the turbines. A summary of the percentage increase in power for all combinations of layouts and terrain can be found in Table 3.
The wind speed contour plots corresponding to the four inflow directions can be seen in Figure 7.   Table 3. The percentage increase gain by using different farm layouts applied to the flat and skewed terrains. All values are normalized by the power output from the gridded layout on the flat terrain. The turbine configuration in the layouts can be found in Figure 6 wake recovery. Two of the corner turbines are also slightly offset to further reduce wake effects from other corner turbines, (e.g. turbines 1, 2, and 4). Turbine 1 and 4 are bound by the site constraint and cannot move further away from the center of the domain. Another interesting feature to point out is that several wakes curve around the hill and other turbines. Turbine 3 shows significant wake curvature when the inflow is from the north and turbine 4's wake bends around the hill. This optimization was performed with both the analytical definition of the skewed Gaussian hill as well as an interpolated version of the hill's mesh coordinates read from file. Both terrain representations resulted in nearly identical optimized layouts.

Conclusion
Complex terrain has a substantial impact on the flowfield within a wind farm, and as our results show, a substantial impact on layout optimization results. We have modified NREL's WindSE RANS code to simulate flow in complex terrain and shown that this leads to significant alterations in the trajectory of wind turbine wakes that are not captured by traditional engineering flow models. Specifically, we observed significant spanwise and vertical velocities at hub height, along with terrain-induced speedups on ridgetops and flow separation on the leeward side of terrain features. These results were validated by comparison against SCADA data from a real wind farm. In our study on 16 turbines, the terrain-aware model enhancements improved the agreement with the normalized power output from SCADA data by 34.8% in terms of mean average error. After validating the model in complex terrain, we performed two layout optimization studies. The first study on a simplified Gaussian hill confirmed that WindSE's automatic differentiation correctly accounts for terrain effects in computing gradients for optimization. In our second  Figure 7. This figure shows the minimal wake interactions for the layout optimized for the skewed hill for all inflow directions. All plots display the streamwise velocity located 80 m above the ground and inflow speed is 8 m/s at hub height on the inflow plane. layout study we considered 9 wind turbines, a skewed hill for terrain, and 4 inflow directions. The optimization strategy discovered in these studies concentrates 5 turbines on the plateau where winds speeds are maximized, while keeping the remaining 4 turbines near the corners. The turbines on the hilltop are placed to take advantage of high winds while avoiding waking each other in the 4 inflow directions we considered. The turbines in the corners maximize spacing in the streamwise direction to allow for wake recovery, but are also slightly offset from each other to avoid a direct wake. Additionally, the optimized layout has many symmetric characteristics that reflect the evenly weighted inflow directions and skew-symmetric terrain. Similar symmetry has been reported in previous WindSE layout optimization studies [1].
These layout optimization results yield a 123.2% improvement in our objective function relative to the starting configuration with a gridded layout. We attribute this large improvement to the fact that the starting configuration only has one turbine on the hilltop where velocities are increased approximately 20%, while the optimized layout places five turbines. The turbines are operating in Region 2 in our model, so the terrain-induced speedup leads to significant power gains. Furthermore, the gridded layout produces substantial wake losses in the four inflow directions we considered, while the optimized layout minimized wake interactions. Additionally  12 our optimization can be performed using interpolated terrain data, which allows for the possibility of optimizing on arbitrary terrain.
These initial simulation and optimization results open many new paths for future research on wind farm performance in complex terrain. We are pursuing further research on optimization of larger wind farms and considering additional inflow directions and speeds to complete a full annual energy production optimization. The impact of complex terrain on the curved trajectory of wind turbine wakes will also likely impact wake steering strategies that are an emerging tool for wind farm controls. Additionally, the vertical motions induced by the terrain will be affected by the thermal stability of the atmospheric boundary layer, making stability an important consideration for further research in complex terrain. Finally, since gradient based optimization methods find local minima but are not guaranteed to reach the global minimum, future studies are needed to assess the impact of multiple restart strategies on the optimization result.