Large-eddy simulation of spectral coherence in a wind turbine wake

This work is mainly dedicated to the study of the characteristics of spectral coherence of turbulence fluctuations in wind turbine wakes. A computational fluid dynamics (CFD) code has been implemented using a large-eddy simulation (LES) approach, which is thought to be conceptually more suitable for studying the turbulence evolution in a wind turbine wake. Comparisons with experimental data from the Nørrekær Enge II Windfarm, in Denmark, and with an analytical model proposed by Panofsky and Dutton have been performed, and the results are found to be in reasonable agreement with both.


Introduction
Due to the cost of land and civil works, wind turbines tend to be located close to each other in wind farms, creating interference effects. Consequently, wind turbine wakes are an interesting topic of study, because the momentum defect and the increased level of turbulence created by turbines in a wind farm may cause a reduction in power output and unsteady loads on other machines. Specifically, turbulence may increase in turbine arrays sufficiently to cause measurable damage due to fatigue and dynamic loads, that have to be taken into account adequately; see Frandsen and Thögersen (1999). Only a direct simulation of turbulence will be able to give us the full knowledge of the turbulence characteristics in the wake. In industrial or environmental applications, where Reynolds numbers are usually very high, direct-numerical simulations (DNSs) of turbulence are generally impossible because the very wide range that exists between the largest and smallest dissipative scales cannot be explicitly simulated, even in the largest and most powerful computers. A large-eddy simulation (LES) approach is conceptually more suitable for studying the turbulence evolution after a wind turbine, because the larger scales of turbulent motion, those that control turbulent diffusion of momentum and energy, are computed explicitly, whereas only the effect of the small-scale motion that tends to be more isotropic and dissipative has to be modelled. Different approaches, reviewed by Crespo et al (1999) and Vermeer et al (2003), have been adopted to study wind turbine wakes. The work presented here is a continuation of the UPMWAKE model proposed by Crespo et al (1985) and Crespo and Hernández (1989), based on the k-ε closure methods, and the explicit algebraic model for the components of the turbulent stress tensor proposed by Gómez-Elvira et al (2005). In all these methods, the Reynolds average over all turbulence scales was used. LES will reproduce the unsteady oscillations of the flow characteristics over all scales larger than the grid size; consequently, a greater detail of the turbulence characteristics is expected to be obtained. An LES model with simplified boundary conditions has been implemented in a CFD code to simulate and characterize the turbulence generated by the presence of a wind turbine. In this work we have performed an LES calculation of the wake by placing the simplified turbine in an environment with turbulence properties similar to the ones of the atmosphere. The wind turbine is modelled by drag forces applied at the grid points corresponding to the projected disc area. This model has been previously applied by Jimenez et al (2007) to study the general characteristics of the Reynolds stress tensor in the wake and make comparisons with experimental data from the Sexbierum wind farm, in The Netherlands, and with analytical correlations. The results were found to be in good agreement with both, suggesting that LES is a potentially useful tool in the investigation of detailed wake flow.
One point of interest that an LES simulation allows one to investigate is the spatial coherence of turbulence variables of the wake. Turbulent velocity fluctuations in different positions in space are correlated, depending on the characteristic scale of turbulence. This correlation will obviously have some influence on the loads affecting structures of some spatial extent, subjected to the turbulent wind. Traditionally, the spectral coherence of the wind speed fluctuations at two points has been used as a quantitative measure of the amount of correlation as a function of frequency (or of spatial scale).
The behaviour of the coherence is quite well known from many field experiments in equilibrium turbulence, but not much is known about coherence in the non-equilibrium turbulence that exists in a wind turbine wake, where the turbulence statistics can be seriously altered by the local generation of turbulence caused by the strong shear introduced by the extraction of energy from the mean flow as shown in Højstrup and Courtney (1993). The turbulence generated by the wake shear will typically appear at length scales of the order of the cross dimensions of the wake, which are much smaller than the typical scales of horizontal turbulence in the atmosphere.
Unfortunately, two different definitions of the coherence exist in the literature. The definition of coherence employed here will be where Q( f ) is the quadrature spectrum, Co( f ) is the cospectrum, and S 1 ( f ) and S 2 ( f ) are the power spectra. The other definition used, sometimes called the 'squared coherence', is given by the square of equation (1); see e.g. Bendat and Piersol (1966). In this work we present results obtained using an LES approach, as well as comparisons with the experimental data given in Højstrup (1999) and with a very simple analytical model proposed by Panofsky and Dutton (1984).

Governing equations
We consider a three-dimensional airflow of incompressible and viscous fluid moving between to parallel planes (channel flow). It is assumed that there is neutral atmospheric stratification and that density variations are negligible. In a DNS, the governing equations consist of the continuity and Navier-Stokes equations, as follows: where the subscripts i , j = 1, 2, 3 correspond to the streamwise (x), vertical (y) and spanwise (z) directions, respectively. In the above equations, u i is the instantaneous velocity component in the i direction, p is the instantaneous pressure, ρ is the density and ν is the kinematic viscosity coefficient.
In an LES simulation, the flow variables are divided into a grid-scale (GS) part and a subgrid-scale (SGS) part by the filtering operation. Then, the filtered continuity and momentum equations are given by whereū i is the instantaneous filtered velocity component in the i direction, andp is the instantaneous filtered pressure. The effect of the unresolved SGSs appears in the SGS stress is as follows: which must be modelled. In this study, τ (SGS) i j is parameterized by the eddy viscosity assumption of Smagorinsky (1963), through the following constitutive relations: where δ i j is the Kronecker delta, ν T is the eddy viscosity, S i j is the resolved strain rate tensor, |S| is the magnitude of the resolved strain rate tensor, is the grid-filter width, which is a characteristic length scale of the largest SGS eddies, and x , y and z are the grid spacings in the three axis directions. In the original Smagorinsky model (Smagorinsky 1963, Lilly 1965, the dimensionless parameter C S is an empirical constant, whereas in the dynamic Smagorinsky model (Germano et al 1990) this parameter is computed dynamically during the solution, making it a function of space and time. The dynamic model is adopted in this work.
Finally, defining a pseudopressure in the following way: the governing equations present the form Temporal averages will be indicated by a bracket, · , and their fluctuations, that is the instantaneous filtered value of a flow magnitude minus its time average value, are indicated by a prime, :φ = φ +φ , where ϕ is any of the flow variables, velocity or pressure. As the average process is steady, time averages will also be the classical Reynolds average. In the following we will always operate with filtered values, and for simplicity we will omit the bar. The subgrid-scale (SGS) stresses (see equation (6)) are found to be negligible compared with the resolved ones | τ (SGS) i j | < 0, 01ρ u ı u j | for the grid used. Thus, it is justified to assume that most of the energy and momentum transport rate is contained in the resolved fluctuations; therefore only the filtered values are considered in the following discussion.

Numerical scheme
A set of channel flow simulations have been performed using the CFD code elaborated at the Center for Turbulence Research, Stanford University by Charles D Pierce (see Pierce andMoin 2001, Wall et al 2002). The numerical method used to solve equations (13) and (14) is a semi-implicit iterative solution procedure that is both economical but also retains most of the stability and accuracy properties of the fully implicit scheme. It is conveniently described in those references. The calculations have been performed with a grid of 256 × 96 × 96 points in the x, y and z directions respectively, which, considering the size of the computational domain described below, ensures a good quality in the shape of the cells. In the vertical direction, the grid points are non-uniformly distributed; in particular, they have been clustered around the turbine. The first grid point is placed at a distance of 0.018h over the ground, where h is the wind turbine hub height. Grid convergence has been verified, by checking that time average results do not change significantly with an increase of the number of grid points.

Flow model
The configuration to be studied is shown in figure 1. The computational domain is a rectangular channel of the shape of a parallelepiped whose sides have lengths L x = 34.9D, L y = 5.6D, and L z = 10.7D. Inside the box there is the wind turbine, which is simulated by several momentum sinks. On the side walls of the channel, z = ±L z /2, periodic boundary conditions are applied. For the upper wall, a boundary condition giving a fixed value of the shear stress is imposed. Periodic boundary conditions are also applied at the entrance plane, x = 0, and the exit x = L x . Thus, the flow field that is represented corresponds to an infinite line of wind turbines separated by the distance L x . It could be argued that if the interaction is really equivalent to that of an infinite number of turbines, the periodic lateral boundary conditions will work only until the wake impinges on the boundary, and consequently this condition will never be satisfied, unless the wake effect disappears in the distance L x . Calculations have been carried out both with and without turbines, and it is found that the length L x is large enough that the flow field incident over the wind turbine is equal to the unperturbed flow without any turbine. It is expected that this incident flow field will have the characteristics of a real atmospheric flow. In particular, that the average flow velocity in the x direction will be given by the classical law of the wall, which in the case of rough wall is expressed as where u * is the friction velocity, κ is the von Karman constant (κ = 0.4) and y 0 represents the terrain roughness. We will check later whether this law of the wall is really satisfied.

Wind turbine model
The next step is to introduce inside the airflow a model of wind turbine that allow us to understand how the presence of the turbine creates a wake and modifies the flow characteristics. The circular disc is approximated by a set of rectangular cells in a Cartesian grid. The model consists of the disc formed by the set of cells (see figure 1) in which we introduce a force in the x direction, proportional to the area of each cell and to the square of the incident average velocity in the central cell: C T represents the thrust coefficient of the wind turbine, A is the frontal area of the cell, and U 0 is the unperturbed velocity of the incident flow in the centre of the disc. The minus sign means that the force has opposite direction to that of the mean flow, and therefore, decelerates the incident stream. Such a model of wind turbine is like a porous disc, which occasionally has been used to simulate wind turbines. Larsen et al (2003) present the results from 3D CFD computations, indicating that the near wake field may be separated into two regions, according to whether or not the flow structures around the turbine blades have developed to an azimuth invariant stage. The extent of the inner region (in which blade flow structures can be identified) is only of the order of one rotor diameter.
The tower shadow could also be simulated with a distribution of sinks; however, as the frontal area of the tower is much smaller than that of the disc, and the expected drag coefficients of the tower and disc are quite similar, the tower effect has been neglected. In the present study, the real blades of the turbine are not taken into account because we are not really interested in the local properties of the flow directly interacting with wind turbine blades, but in the downstream evolution of the flow characteristics. In particular, we principally want to study the turbulence added in the wake. That turbulence is mainly generated in the shear layer of the near wake, as is shown in figure 2, where a schematic view of the horizontal cross section is shown. As can be seen in this figure, there is a core region where the flow is decelerated by the wind turbine, and outside it the velocity has the ambient value; the difference between these two velocities creates the shear layer that has a ring-like shape. Turbulent diffusion makes the shear layer thickness increase with downstream distance, and, at a certain distance downstream (about two to five diameters), the shear layer reaches the wake axis. It should be taken into account that this shear layer in the near wake is made up of concentrated vortices, and that the flow is not really turbulent outside them. The superficial resemblance to classical turbulent shear layers, as in wakes and mixing layers, occurs largely because wind turbine wakes are mostly viewed in an average sense: the flow at any point is treated as an ensemble average over many blade revolutions, so that azimuthal variations in the mean velocities, as seen by an observer rotating with the blades, appear as 'turbulence' to a non-rotating observer. It could be thought that, in addition to the linear momentum extraction in the actuator disc, a rotational speed must be added in the wake to compensate for the torque made on the turbine, in order to achieve the angular momentum conservation. Surely this rotational component must exist, but it can be shown that it is very small compared not only with the mean velocity of the flow, but even compared with a characteristic turbulent velocity. From the actuator disc theory, an order of magnitude calculation based on angular momentum conservation yields ρ π D 2 where a is the induced velocity factor of the actuator disc, and M the torque generated in the machine. In the case of the wind turbines Nordtank NKT with 300 kW of rated power for a wind speed of 15.5 m s −1 , and a diameter of 28 m, the angular speed of the rotor is about 7.5 rad s −1 , and then the torque can be set approximately to 40 000 N m. Using the value a = 1/3 we get v θ ≈ 0.35 m s −1 . That is of the order of 2% of U 0 , being the turbulent oscillations of velocity components in the range of 10-15% of U 0 .

Boundary conditions
The boundary condition in the bottom side of the computational domain, y = 0, (see figure 1) is expressed by the following equations, frequently used in meteorology: and are proposed by Senocak et al (2004). Very similar approaches are used by Moeng (1984) and by Piomelli and Balaras (2002).
In expression (19), τ iy are the instantaneous global shear stresses and d is the distance of the first node from the wall. Note that condition (19) is consistent with the law of the wall, equation (15), although in that equation the temporal average values of velocity and shear stress were used, whereas in the previous condition we are operating with instantaneous values. In particular, τ xy (x, y) = ρu * 2 is the classical definition of u * in the law of the wall.
In the top side of the domain, y = L y , we adopt Equation (23) guarantees a uniform value (that means no y-variation) of the average shear stress in the regions that are not affected by the wake.

Experimental data
In Højstrup (1999) there is presented a set of results related to spatial coherence of turbulence obtained from measurements done in a wind turbine array, the Nørrekaer Enge II Windfarm. It contains 42 × 300 kW Nordtank turbines which have a rotor diameter of D = 28 m and a hub height of h = 31 m, and it is located in North Jutland, Denmark. The configuration of the park consists of two sets of three rows of turbines with seven turbines in each row; see figure 3(a). There are two masts in the wind farm. One of them, called M2, is sited deep inside the park. It is placed on the line between turbines E5 and F6 at a distance of 55 m from F6. This mast experiences wakes from different turbines. When the wind direction is about 30 • , the mast is inside the wake generated by F6 at a distance of 2D from it, and when the wind direction is the opposite, the mast experiences the wake created by E5 at a distance of 7.5D. Two cup anemometers were mounted on the opposite sides of this mast ( figure 3(b)), such that coherence measurements could be made at 2D downstream from turbine F6 and 7.5D downstream from turbine E5. These experimental data were subdivided into wind speed categories, 6-8 m s −1 , 8-10 m s −1 and 10-13 m s −1 , called respectively '7 m s −1 ', '9 m s −1 ' and '11.5 m s −1 ' cases. The wind speed used for the selection was the wind speed in the wake measured at the highest position, which approximately corresponds to h + D/2.

Panofsky and Dutton's model
Experimental and numerical results for the lateral and vertical coherence of the longitudinal velocity component will be compared with a simple model given by Panofsky and Dutton (1984) and also applied by Højstrup (1999): where f is the frequency, s is the separation, U the mean speed and a = 6 + 11(z 2 − z 1 )/(z 2 + z 1 ).

Results from LES calculation
A set of LES calculations has been performed to simulate the conditions of the experiments described above, and coherence was calculated at points sited in the wake at 2D and 7.5D distance downstream from the turbine. Their relative positions are the same as shown in figure 3(b). Also, the lateral and vertical coherence of turbulence in the absence of wake was obtained from simulation. The parameters needed by the model were set as follows. In Højstrup (1999) it is specified that the terrain in which the park is placed is an old seabed and extremely flat, about 1 m above sea level and with expected roughness lengths of a few centimetres. Then, for the ground roughness we adopted y 0 = 3 cm. The thrust coefficients C T are obtained from the thrust curve of the Nordtank 300 kW turbine, leading to the values 0.68, 0.63 and 0.56 for C T when incident velocities at height h + D/2 are respectively 7 m s −1 , 9 m s −1 and 11.5 m s −1 . For all the three cases simulations lasted 8192 s, which permits us to obtain eight series of 1024 s each. After FFT computation was done, cross-spectra were calculated and band averaged (bandwidth 0.07 decades), upon which crossspectra in each category (wind speed and height over the ground) were averaged together. Then, the coherence was evaluated from equation (1). All the procedure is the same as that used by Højstrup (1999) when obtaining the coherence from experimental data.

Comparison of the characteristics of the incident
unperturbed flow over the turbine with those of a typical atmospheric flow. In figure 4(a) there is given the calculated velocity profile incident over the turbine, and it is compared with the logarithmic law expressed in equation (15); the agreement is good in the zone where the wind turbine is placed. In the upper part of the domain the calculated value of the velocity is larger than the one given by the law of the wall. This is due to the boundary conditions (22)-(24), imposed in the upper side of the computational domain. With other boundary conditions the velocity profile was closer to the logarithmic one, but the components of the turbulent stress tensor were poorly reproduced. In future works we will try to correct this effect. In figure 4(b) there are shown the distributions of all the calculated components of the Reynolds stress tensor τ i j = −ρ u i u j . For the grid used in calculations, the SGS stresses τ (SGS) i j are found to be completely negligible compared with the resolved ones, except for the component τ xy in a very small region close to the walls, since at them u y = 0 at any time, and then u y = 0. Nevertheless, this does not happen in the vicinity of the actuator disc because it is not treated as a solid disc, but as a porous one; then velocity fluctuations are not null close to the simulated rotor. The components τ xz and τ yz are zero because of symmetry. As indicated in equation (21), the mean turbulent shear stress of the incident flow should be τ xy = ρu * 2 , which is in excellent agreement with the calculated value. The three diagonal components are somewhat smaller than the values estimated by Panofsky and Dutton (1984), from a set of experimental results, which are also indicated in international standards: However, the numerically predicted values of the longitudinal component of turbulence are not too far from those measured in the Nørrekaer Enge II Windfarm, as can be seen in figure 4(b). The coherence of the unperturbed flow will be examined in the next section.    • On the other hand, both the numerical and the experimental results predict values of the vertical coherence that are larger than those given by the model for large frequencies.   Figure 11. Vertical coherence for 13 m separation and U = 7 m s −1 in no wake, 2D wake and 7.5 wake cases. The straight line represents the model given by equation (25).  Figure 13. Vertical coherence for 13 m separation and U = 11.5 m s −1 in no wake, 2D wake and 7.5 wake cases. The straight line represents the model given by equation (25).
same band width as Højstrup (1999) with the purpose of making a proper comparison. Alternatively, these tails may also be explained because there are systematic errors on coherence estimates when the number of series (or realizations) is limited; in our case this number is eight (of 1024 s each). Højstrup's data also present a limited number of realizations. From Mann (1994) (equation 4.5), it can be seen that if the real coherence is 0, then the one estimated with only eight realizations is rather significant (more than 0.3). • In both experiments and calculations, there are not big differences between the vertical coherence inside and outside the wake, except for figure 12, U = 9 m s −1 , where the experiments show an increase of coherence in the near wake (2D). • On the other hand, the lateral coherence presents smaller values in the near wake, this effect being more noticeable in experiments. See figures 5-10. • No significant effect of the wake is noticed in the far wake (7.5D), where the coherence seems to be similar to the upstream values (no wake).

Conclusions
A simplified LES model is proposed to simulate the turbulent flow in the wake of a wind turbine. The turbine is simulated by a set of local sinks of momentum distributed across the rotor disc, without reproducing the blade details. The turbulence characteristics, in particular the coherence at every pair of points of the computational domain, can be obtained. Comparison has been made with values of vertical and lateral coherence measured by Højstrup (1999) and an analytical correlation established in the literature for free atmospheric turbulence. A reasonable agreement is found. The interpretation of the tendencies of the coherence behaviour is not straightforward, and it will be further elaborated in future work.
These results indicate that this LES model, with the simplified approach to simulate the rotor, is a very useful tool to simulate real turbulent characteristics in wakes. Other turbulence characteristics such as spectra, limit values of flow