Theory and Verification of a new 3D RANS Wake Model

This paper details the background to the WakeBlaster model: a purpose built, parabolic three-dimensional RANS solver, developed by ProPlanEn. WakeBlaster is a field model, rather than a single turbine model; it therefore eliminates the need for an empirical wake superposition model. It belongs to a class of very fast (a few core seconds, per flow case) midfidelity models, which are designed for industrial application in wind farm design, operation and control. The domain is a three-dimensional structured grid, a node spacing of a tenth of a rotor diameter, by default. WakeBlaster 5 uses eddy viscosity turbulence closure, which is parameterized by the local shear, time-lagged turbulence development, and stability corrections for ambient shear and turbulence decay. The model prescribes a profile at the end of the near-wake, and the spatial variation of ambient flow, by using output from an external flow model.


Theoretical Background
The WakeBlaster wind farm simulator is based on a Reynolds-Averaged Navier-Stokes (RANS) set of equations, which is used to solve the propagation of wake dissipation through the farm domain, in Cartesian 3D coordinates. In order to account for the fluctuation term of the velocity vector, it uses eddy viscosity turbulence closure, where the eddy viscosity is calculated from the combined wake and ambient wind speed shear profiles.

RANS Equations
The wake model uses RANS equations for momentum conservation, and mass flow conservation to calculate the three components of wind velocity in the axial, lateral and vertical directions. Cartesian 3D vectors are used for displacement − → x and wind speed relative to ambient − → u : − → x = x, y, z − → u = u, v, w , where the first element of the vectors (x) is the streamwise component, the second element (y) is horizontal and perpendicular to (x), and the third element (z) is vertical (starting from 70 the ground up) and makes up a right-hand coordinate system.
The Reynolds averaged momentum and mass conservation equation can be expressed in two dimensions, for either a free jet or a wake submerged in an incompressible fluid, as given by Abramovich (1963): representing the momentum in the flow direction, where u , v and w denote fluctuations from mean values, and ν the viscosity 75 and ρ the density of the fluid. The corresponding continuity equation is The momentum equations in transversal directions are not considered in the description of a free jet or wake present beyond the near wake of a wind turbine.

80
The following simplifying assumptions are applied by Abramovich for a stationary free wake, expanding into an infinite region: Viscosity The effect of molecular viscosity is small ν ∂ 2 u ∂y 2 = 0 compared to the turbulent viscosity Pressure Flow pressure gradients can be neglected in most cases 1 ρ ∂p ∂x = 0 Stationary The flow is stationary with respect to the mean velocities ∂u ∂t = 0 Thin shear layer approximation Fluctuations along the flow change much slower than in the transversal direction ∂u u ∂x = 0

85
After substituting the continuity equation and applying the simplifying assumptions Abramovich (1963) obtains: or expanded to three dimensions: Using the Boussinesq eddy viscosity assumption, the stress components u v and u w are expressed as: where denotes the isotropic eddy viscosity. The streamwise variation in transversal velocities ( ∂v ∂x and ∂w ∂x ) is small compared to the transversal variation of streamwise velocity ( ∂u ∂y and ∂u ∂z ). The spatial variation in eddy viscosity can be neglected in first approximation and is therefore approximated as a constant, the governing momentum conservation equation can now be written as:

Numerical Solution
The ambient wind field is determined by an external flow model, and it determines the inflow conditions and spatial variations over a site. The turbine is represented by its hub height, diameter and other readily available and measured characteristics.

100
The waked wind field is set up by creating a two-dimensional flow plane, which forms a cross-section along the y and z axes of the velocity vector u. The flow plane is propagated downstream along the x coordinate and when it passes a turbine, a wake is injected into the flow plane.
The grid spacing is set by default to a tenth of rotor diameter. In the vertical direction, the grid starts at the ground z = 0 and reaches up to a default height of three rotor diameters or 31 grid layers. In the horizontal direction, the grid is expanded, as 105 required, to enclose each wake injected into the flow plane with an additional four rotor diameters to the side to allow for wake expansion.

Wind Turbine Momentum Extraction
Axial-momentum theory prescribes pressure building up in the induction zone upstream of any wind turbine or wind farm, and pressure recovery in the near-wake downstream of the rotor. The momentum that each of the turbines extracts in the process is the wind speed dependent thrust coefficient, as a function of the idealised incident wind speed, U inc , at each turbine location, without the presence of the turbine.
In the model, the momentum deficit is injected at the end of the near-wake (which is assumed to be at 2 diameters downstream of the rotor) of each turbine, and it is distributed over an expanded rotor area, using the blunt bell-shaped wind speed deficit profile from Lissaman et al. (1982). The centre-line wind deficit relative to incident wind speed D m , experimentally determined 115 by Ainslie (1988) at a downstream distance of 2 diameters is used as a function of inflow turbulence I inc and thrust coefficient The radial width of the profile is then derived by ensuring momentum conservation with regard to the thrust coefficient of the turbine. 120

Flow Plane Propagation
The flow plane is propagated according to equation 6 using the alternating direction implicit (ADI) method described by Peaceman and Rachford (1955);von Rosenberg (1983), where it is alternately solved in the xy and xz planes, incrementing the x (downstream) coordinate by half a propagation step between each solving plane, so that both planes are solved once per step. By solving for each row or column in the flow plane, and by employing the central difference method, the problem 125 can be set up numerically in a tridiagonal matrix equation, which can then be solved efficiently for the axial velocity, u, by the Thomas algorithm Thomas (1949), described for example in Burden and Faires (2001). In 3D Cartesian coordinates the In practice, due to the assumption of incompressibility, this formulation will lead to a local velocity shear, resulting in nonzero lateral and vertical velocities that are infinitely far from the source of shear. In reality this would not be the case, due to the compressibility of air. Therefore, in order to account for the effect of compressibility, a spatial damping term is introduced, so that v and w tend to zero at y = −∞, y = ∞ and z = ∞: where γ is a user-configurable positive constant that determines the strength of lateral and vertical velocity damping. As these integrals are indefinite, boundary conditions must be assigned. In the vertical direction, it is given that vertical velocity at ground level is zero, as mass flow cannot pass into or out of the ground. Therefore, the condition w z=0 = 0 is applied, leading to: In the lateral direction, the physical boundary conditions are that v y=−∞ = v y=∞ = 0, because the wind farm wakes cannot induce lateral velocity far from the farm. However, for numerical purposes, the size of the flow plane is constrained, and it cannot be guaranteed that the velocity will reach zero on both sides of the flow plane. Therefore, the lateral velocity is integrated in each direction, starting from zero, and the mean of the two is taken. This is expressed as: where y min and y max are the lateral location of the edge of the flow plane.

Eddy Viscosity Calculation
The key term controlling the rate of wake dissipation is eddy viscosity. Eddy viscosity has dimensions of length squared over time, and it can be represented by multiplying a length scale of the shear layer by a velocity scale of the flow field.
WakeBlaster calculates eddy viscosity from the shear profile of axial velocity in the yz plane. In order to do this, it creates 155 a combined flow plane of the ambient wind speed, U amb , multiplied by the solved wake flow plane, u, which is relative to ambient wind speed. In neutral atmospheric conditions, the ambient wind speed is calculated as a logarithmic function of height above ground: where u * is the friction velocity, taken to be 2.5 times the value of standard deviation of the axial wind velocity, κ is the 160 von-Karman constant (value = 0.4) and z 0 is the roughness length. The unknown parameters are determined from inputs to the simulation, such as wind speed and turbulence intensity at a particular height (usually the hub height of one of the turbines).
The eddy viscosity is then calculated for every point in the flow plane, using the following process: 1. Create a combined flow plane by multiplying the ambient surface layer wind speed profile by the waked flow plane velocity u. 3. In each of the two directions, the component of eddy viscosity is calculated as i = ∆u i Λ i , where ∆u i is the difference between minimum and maximum velocity and Λ i is the distance between the maximum and minimum points. This 170 process is shown in figure 1.
4. The overall eddy viscosity is the calculated as = k 2 y + 2 z , where k is a positive calibration constant, which although configurable is considered to be independent of wind farm size and layout. For a logarithmic wind speed profile in the vertical direction with no lateral variation, this method leads to an eddy viscosity that is proportional to the height above ground.

Eddy Viscosity Lag
The eddy viscosity, as so far described in section 2.4, is solely based on the wind shear profile. However, no newly created shear profile instantly generates turbulence, and therefore eddy viscosity -in reality, there is a lag between the change in a shear profile and its effect upon eddy viscosity and wake dissipation. In WakeBlaster this lag is formulated in terms of downstream distance, and it has two distinct models.

180
The 'fixed' model obeys a first order lag equation: where is the lagged eddy viscosity, Λ is the length-scale defined in the previous section, and a positive constant defining the lag length relative to the length scale and is considered to be independent of wind farm size or layout.
The 'turbulence dependent' model gives a larger lag distance when the eddy viscosity and turbulence are low, and it obeys 185 the following equation: where φ is a positive parameter that determines the strength of turbulence on the lag length, and λ max is also a positive parameter that corresponds to the lag length when turbulence is zero. Both parameters are calibrated against extensive wind farm observational data and are considered to be independent of wind farm size and layout.

Atmospheric Stability
When simulating atmospheric conditions that are not neutral, the calculation of eddy viscosity is modified. This modification uses the Monin-Obukhov length, L, and the concept of non-dimensional wind shear, φ m , which is defined by Businger (1971), as: Furthermore, according to Businger (1966), the non-dimensional wind shear is empirically approximated as what tends to be known as the Businger-Dyer relationship: where ζ = z L . The ambient wind speed shear profile is then modified by introducing ψ m : 200 where: where ζ 0 = z0 L . Furthermore, the vertical component of the eddy viscosity, z , is also modified by the non-dimensional wind shear: The horizontal component of eddy viscosity y , is left unmodified.

Wind Turbine Power Calculation
WakeBlaster calculates the power output using power curve input from the user. In order to calculate accurate power, corresponding to the variant wind speed across the rotor, a rotor equivalent wind speed (U rot ) is calculated. This is done by first calculating the combined ambient and wake axial velocity (U = U amb u at the rotor plane), and then integrating across the rotor 210 disk area: where n is an integer. A popular approach is to use n = 3 as suggested in IEC61400-12-1 (2017), based on the principle that the power available in the wind is proportional to the cube of the wind speed. However, WakeBlaster uses a value of n = 1 by default as turbines will not be able to realise the full potential of a sheared inflow over the rotor. As this method is performed on 215 the combined ambient and wake axial velocity, the effects of wind shear on power production are implicitly included whenever the severity of the wind shear depends on the turbulence and atmospheric stability of the flow case.
A general directional variability of the wind within each flow case is included in a standard power curve. A rotor yaw angle can be set per turbine, to consider in the power calculation a known average directional misalignment with the rotor plane. A model to modify the power curve for site specific directional variability over the rotor, for example changes with height or for 220 specific meteorological conditions, is not included in the model.
WakeBlaster uses IEC methods in IEC61400-12-1 (2017) to adjust the power curve for air density and turbulence intensity.
The rotor equivalent turbulence intensity is also calculated using the integral method above, but instead using a value of n = 2.

Verification
In this section, the grid dependence and sensitivity is analysed and an estimate of the numerical uncertainty is thereby pro-225 vided. Computational performance for large wind farms is verified, and offshore wind farm model predictions are inspected graphically, for plausibility.

Grid Dependence and Sensitivity
The model uses a structured grid, in terrain following coordinates. The grid resolution is scaled with a length scale characterising the specific flow -the rotor diameter. The grid is equally spaced in all directions, and no stretching, compression, or nesting 230 is applied to any part of the domain. The minimalist design is computationally efficient, and it avoids potential numerical errors -at grid interfaces which do not match, for example.
The solver is designed for a single purpose: to model the impact of wind turbines on the underlying flow and the consequential wind farm wake losses. A wind turbine's wake scales with its rotor diameter and its height above ground. In order to match Analysis of the sensitivity of model results to changes in grid resolution verifies that the results are not sensitive to grid resolution over the expected range of application. Challenges could arise -for example, when using an average resolution in wind farms with mixed turbine diameters and turbines mounted at low hub heights. In an annual energy calculation, the overall wake loss is composed of several thousand individual flow cases. Wake loss model errors are commonly estimated to be in the 240 range of 10-20 percent, relative to the average annual wake loss. Numerical errors should be one order of magnitude lower.
Ignoring error compensation between flow cases, an error of 1-2 percent (relative to the wind speed difference for an individual flow case) is acceptable.
The grid dependency study was carried out for the following scenario: A single turbine, with ambient wind speed perpendicular to the rotor plane. Ambient conditions were a wind speed of 8 m/s, neutral atmospheric stability, and a turbulence 245 intensity of 10%. The wind turbine type (V100-1.8) is described by its geometry and thrust coefficient, which is c t = 0.8 for the scenario. The key target value investigated is the wind speed relative to ambient wind speed, at hub height and at several distances downstream.
The sensitivity was tested in a flow case with a strong wake, and the results are presented in 2. The error in wind speed is presented relative to a hypothetical error value, which is calculated using a Richardson extrapolation for an infinitesimal grid 250 spacing, as suggested by Roy (2003) for mixed order numerical schemes.
The numerical error, due to grid spacing for an operational range of up to just above 0.1 diameters, is below 1 %. At a coarser resolution the model can no longer resolve the structure of the flow sufficiently. The current choice of grid resolution (0.1 diameters) represents a reasonable compromise between computational efficiency and model accuracy.
The grid resolution in the model scales automatically with the rotor diameter. Neither the grid nor the resolution are variables 255 which should (under normal circumstances) be adjusted by any user.

Computational Performance
WakeBlaster is a medium-fidelity tool, which is typically capable of running each flow case in a few seconds, on the single core of a modern processor. With the default settings (a flow plane resolution of 0.1 rotor diameters and a domain height of three diameters), the time (in seconds) to run a single flow case (T f c ) is (on an Intel i5 8 th generation processor) approximately: where A is the area of the wind farm and D is the rotor diameter. The T f c is proportional to the area of the wind farm, and (at equal turbine density) to the number of wind turbines in the wind farm. However, the exact time will depend on the wind farm's layout, the wind direction, and the architecture of the processor. The T f c is also proportional to the cube of the flow plane resolution, although results do not show any significant improvement in accuracy when the resolution is increased.

265
For example, a typical flow case for Horns Rev -a wind farm with 80 turbines arranged in a grid, with inter-turbine spacing of 7 diameters -runs in about 5 s. Unless hysteresis effects are included in a time series simulation, each required flow case remains independent from the others, allowing many flow cases to be run in parallel. As WakeBlaster is hosted on the cloud, this allows a high level of parallelisation across tens or hundreds of processors, meaning that an energy assessment consisting of (for example) 2,000 flow cases can be completed in a matter of minutes, even for large wind farms.

Visual Verification
Using a three-dimensional wake model, it is possible to create plots of the three-dimensional waked flow field for the complete wind farm, for a particular flow case. This article presents a visualisation of a single flow case from the Lillgrund wind farm, located in the Øresund Strait, between Sweden and Denmark. The Lillgrund wind farm presents a good case study, because the small spacing between turbines (3.3 and 4.3 rotor diameters, along the two principal rows) leads to large wake effects. The 275 layout is shown in figure 3. The turbines have a rotor diameter of 93 m and a hub height of 68 m above mean sea level.
Three cross-sectional slices in the xy, xz and yz planes, for a flow case of 8 m/s wind speed, 270 deg wind direction, 6 % turbulence intensity, and neutral atmospheric conditions, are presented in figure 4.
These simulations indicate that there is significant interaction between wakes originating from individual turbines, and this supports the assumption that the wakes cannot be modelled independently. The wake interaction leads to a complex wake 280 shape downstream of the wind farm. The low hub height of the wind turbines (68 m), relative to their rotor diameter (93 m), results in significant ground-wake interaction effects. As ambient mixing from below is limited, single turbine wakes become asymmetrical in shape, and the point of greatest deficit drifts downwards to below hub height.

Limitations
The code is a mid-fidelity code designed to be fast and capable of simulating projects with several thousand turbines, working 285 with limited amount of readily available input data and be used in an iterative design process. This limits the level of detail that can be included in the sub-models.
-No direct interaction between the turbines and no description of the axial pressure gradient are included in the model.
The induction zones directly upstream and downstream (near wake) of turbines can overlap and interact. This may lead to changes in turbine performance and turbine characteristics and no attempt has been made to quantify such effects.

290
-A basic representation of the the ambient flow is used as input to the model. The wake is modelled as a perturbation of the underlying flow. No attempt has been made to model a two-way interaction with the atmospheric boundary layer.
-The model uses the directional speedups predicted by a suitable flow model (for example in a RSF/WRG format) to account for spatial variation of the wind resource, for example due to orography, or roughness. Further complex terrain effects, like flow separation, are not considered.
295 Figure 4. Plots of the axial velocity in the wind farm relative to ambient wind speed for a flow case of 8 m/s with the wind from due West.
From top to bottom: xy (birds-eye) slice at hub height; xz (side-on); yz (front-on). The white lines show the corresponding planes of the other plots. The xy plot is taken at the turbine hub-height above sea-level -68 m.
terrain or due to meteorological factors), downstream wake locations may not be accurate.
The WakeBlaster model undergoes continuous, data driven improvement, and refined models will be added successively.

Conclusions
This is the first publication to present the theoretical background of WakeBlaster in some detail. WakeBlaster is a recently 300 developed 3D-RANS solver that is specialised to simulate the waked flow field in wind farms. The characteristics of this model show the desired performance balance between speed and level of detail.
Code availability. WakeBlaster calculations are provided as a cloud service and designed for integration in other software packages. Wake-Blaster is available from ProPlanEn directly (www.wakeblaster.net) and through third party implementations. WakeBlaster has been integrated in EMD's WindPro software and is scheduled for release in May 2020.