Reduced order model of the inherent turbulence of wind turbine wakes inside an infinitely long row of turbines

The turbulence in the interior of an idealised wind farm is simulated using Large Eddy Simulation and the Actuator Line technique implemented in the Navier-Stokes equations. The simulation is carried out for an 'infinitely' long row of turbines simulated by applying cyclic boundary conditions at the inlet and outlet. The simulations investigate the turbulence inherent to the wind turbines as no ambient turbulence or shear is added to this idealised case. A Reduced Order Model for the highly turbulent flow deep inside a wind farm is proposed based on a Proper Orthogonal Decomposition. The reconstructed flow is shown to capture the large scale motions of the highly turbulent flow.


Introduction
In large wind farms, both onshore and offshore, the wind turbines in the interior will operate in the wake of the upstream turbines irrespective of the wind direction. A wind turbine operating in the wake of one or more wind turbines yields less power and is subject to greater loadings from the induced turbulence, leading to increased fatique. Understanding this highly dynamic wake interaction is important for wind farm optimization.
Several analytical wake models exist, e.g. Jensen [8] and Frandsen et al. [4]. Generally, these models are based on simple single wake calculations and different assumptions are used to superpose merging wakes to describe the overall wake interaction inside wind farms. These methods are based on steady state considerations, and are developed with the aim of predicting mean quantities, e.g. the mean velocity deficit in the wake. Therefore, the models exclude important details of the dynamic wake interaction and turbulence properties needed to assess the turbine performance and loadings. Frandsen [5] proposed a model for fatique loading on rotors using the effective turbulence intensity as the only governing parameter, arguing that the various changes in turbulence properties usually correlates with the standard deviation of wind speed fluctuations. Here, the effective turbulence intensity is taken as a design variable, and not as a physical quantity. The model has shown improved accuracy compared to the existing engineering models, albeit still occassionally resulting in large discrepancies compared to measurements. The analytical models have also been combined with models describing the so-called meandering or large scale motions of wake. The Dynamic Wake Meandering model assumes that the large scale motions of the wake are directly related to the large scales in the atmospheric turbulence, see Larsen et al. [10] and [11].
Measuring the wake interaction is not a trivial task and measurements contain large uncertainties due to the complex and everchanging inflow conditions. Therefore, numercial simulations are employed, where all relevant properties are described.
The wake interaction in wind farms is examined using Computational Fluid Dynamics (CDF) and Large Eddy Simulation (LES) as well as the Actuator Line method. These methods couple the full Navier-Stokes equations with a blade element approach to include the effect of the wind turbines on the flow without having to fully resolve the boundary layers on the rotor blades. The Navier-Stokes based actuator disc method was introduced by Sørensen and Myken [16] to simulate axisymmetric flows around wind turbines. An extension to full 3-dimensional flows was later introduced as the Actuator Line technique (AL) developed by Sørensen and Shen [17]. These methods have been used extensively for wake studies on single and multiple wind turbines in wind farms, e.g. Troldborg [19] and Ivanell [7]. Similarly, Lu and Porté-Agel [12] and Calaf et al. [3] studied the turbulent flow inside large wind farms in the atmospheric boundary layer using LES. The simulations have provided valuable insights into wake dynamics, but run times are long and thus unsuitable for practical use, e.g. to optimise wind farm layout which requires numerous runs.
The tradeoff between capturing the complex dynamics of the resolved flow and computational speed is still a pending issue when applying wake models to wind farm optimization. The present work includes the analysis of numerical simulations of the complex wake interaction deep inside an idealised wind farms and suggests the foundation for a more dynamic and physically correct wake model. The proposed wake model is derived from a Proper Orthogonal Decomposition(POD), which is an optimal decomposition of the turbulent flow into large coherent structures and optimal in terms of energy content for each mode as opposed to e.g. a Fourier decomposition.
2. Methods 2.1. Numerical Simulation 2.1.1. Flow Solver The numerical simulation is performed using the 3D flow solver EllipSys3D, developed as a collaboration between DTU(Michelsen [14]) and Risø(Sørensen [18]). This code solves the discretised incompressible Navier-Stokes equations in general curvilinear coordinates using a block structured finite volume approach. EllipSys3D is formulated in primitive variables (pressure-velocity) in a non-staggered grid arrangement. In the present work, the pressure correction equation is solved using the PISO algorithm and pressure decoupling is avoided using the Rhie-Chow interpolation technique. The convective terms are discretised using a hybrid scheme combining the third order accurate QUICK scheme and the fourth order CDS scheme. This technique was employed as a compromise between avoiding unphysical numerical wiggles, occurring when using a pure fourth order scheme, and limiting numerical diffusion due to the upwind biasing nature of the QUICK scheme.
The influence of the wind turbine is simulated using the AL technique. The AL technique applies body forces distributed along rotating lines representing the blades of the wind turbine. The actual body forces used in the AL method are calculated using Flex5, a full aeroelastic code for calculating deflections and loads on wind turbines, see e.g. Øye [20]. The aerodynamic loads are calculated using a Blade Element Momentum (BEM) method modified to include the dynamics of the induced velocities of the wake, see Sørensen and Shen [17] for full details on the AL technique. The flow field is thus determined by solving the 3D incompressible Navier-Stokes equations using Large Eddy Simulation and applying the body forces, f . The sub-grid scale (SGS) viscosity used in the LES is modelled using the vorticity based mixing scale by Ta Phuoc et al. [15]. The advantage of representing the individual blades by line-distributed forces is that much fewer grid points are needed to capture the influence of the blades, as compared to simulating the actual geometry of the blades. Therefore, the actuator line model allows for a detailed study of the dynamics of the different wake structures, such as the tip and root vortices, using a reasonably low number of grid points. The drawback of the method is that it relies on the quality of tabulated airfoil data. Airfoil data corresponding to the NM80 turbine is used in the present work. The 2D airfoil data is corrected to account for 3D effects. The NM80 is a three bladed wind turbine with a radius of R = 40m and rated to 2.75MW at V hub = 14m/s. A controller is implemented for the turbine which is a combination of a variable speed Pcontroller for small wind speeds(V hub < 14m/s) and a PI-pitch angle controller for higher wind speeds. Details of how a similar controller operates can be found in Hansen et al. [6]. The implemented controller essentially means that the rotor is not constantly loaded, and as such it operates like a real turbine.

Domain Setup
The wake interaction is examined in an idealised case of an infinitely long row of turbines with initial uniform inflow of U 0 = 15m/s, which is above rated speed for the wind turbine. Uniform inflow means that no atmospheric boundary layer and no ambient turbulence is added, making it possible to study the turbulence inherent to the wind turbines. Furthermore, the situation is idealised by excluding ground and gravitational effects. The 'infinitely' long row of turbines comes from applying cyclic boundary conditions in the flow direction, i.e. coupling the inflow condition to the outflow. Applying cyclic boundaries effectively reduces the domain to contain only three turbines, which continously introduce more and more turbulence into the flow. The simulations have been run using three turbines instead of only one to minimize the effects of the streamwise periodicity of the domain, which potentially could introduce large scale motions.
The cyclic domain consists of a 162 block structure in a Cartesian grid with a total of approximately 18 · 10 6 grid points. The domain is 20R × 20R × 36R in lateral(Y ), vertical(Z) and streamwise direction(X), respectively. The size of the domain is choosen to avoid numerical blockage by adhering to a blockage ratio of 0.8% < 3% as recommended by Baetke et al. [1].
Far field boundary conditions (U = constant) are applied on the transverse and vertical boundaries, i.e. making it a single row of wind turbines, not an infinite wind farm in the transverse direction. The grid is equidistant in the streamwise direction, which is necessary to maintain and model the turbulence correctly as it is transferred using the cyclic boundary conditions. The grid is non-equidistant in the transverse and vertical with a higher resolution in an equidistant region for ±3R. The three rotors are located in (4R, 0R, 0R), (16R, 0R, 0R), (28R, 0R, 0R). Hence, the spatial distance between two consecutive turbines is 12R. Each rotor blade is resolved by 20 grid points.
The nature of the simulations with cyclic boundary conditions and long run times would usually lead to a loss in mass flux as the turbine extracts energy from the flow. However, the mass flux is continously corrected to account for the energy extracted by the turbine, i.e. equivalent to adding a pressure gradient to account for the mass flux loss. Thereby, ensuring that the mass flux within the domain is kept constant at the initial level. This correction is not constant, as the turbine does not operate with a constant thrust due to the controller.

Proper Orthogonal Decomposition and Reduced Order Model
Proper Orthogonal Decomposition (POD) is applied to all three velocity components extracted 1R upstream of the first turbine, i.e. the incoming flow as seen by the turbine. POD is a statistical method yielding an optimal linear subspace, optimal in terms of the variance of the energy content, i.e. it ensures that the included modes contain the most energitic states. POD is a proven method for examining coherent structures for stationary random processes, see Berkooz et al. [2] for a general overview of POD. As such, POD provides a more efficient decomposition of turbulent flows with large scale structures compared to Fourier decomposition. A brief overview of POD as adapted from Jørgensen et al. [9] is given in the following. The velocity components (u = (u, v, w)) of each of the N extracted slices are organised in vectors u i , which constitute a matrix defined as: It is common to subtract the mean value u 0 , hence reducing the number of dimensions by one. Therefore, u ′ j is defined as: The N × N auto-covariance matrix is then given as R = U T U. Defining an eigenvalue matrix with and matrix of orthonormal eigenvectors: yields a eigenvalue problem of the form: The basis vectors or POD modes can be written as: These are spatial modes, which reveal the turbulent structures. Hence, the original vectors u ′ j containing the velocity components of slice j can be reconstructed from where a kj are the elements of A = Φ T U.
Since the eigenvalues or POD modes are sorted, it is possible to approximate the flowfield by truncating the reconstruction by only including the first K POD modes, i.e. where u is the POD truncated reconstruction of u. POD can be considered as an energy filter and as such it is possible to construct a Reduced Order Model ( u) based on the above decomposition. A contains the temporal development of the POD modes, and since the flow has approached a random stationary process the idea is to truncate the reconstruction even further, because the eigenfunctions are composed of Fourier modes. The dominant frequencies of the temporal eigenfunctions can be identified through spectral analysis and by choosing an optimal combination of these dominant frequencies f k,m it is possible to reconstruct the flow field The coefficients β k,m are determined by optimising the reconstruction in a linear least square sense. Here, k denote the POD mode and m denote the frequency number for a given POD mode, sorted in terms of influence after the linear least squares fit.

Results
A total of T = 4292s ≈ 70min has been simulated, which corresponds to passing a total of approximately 130 turbines based on the free stream velocity. Figure 1 shows the spatially averaged (within ±2R) streamwise velocity (and the filtered mean) as function of time for the first of the three turbines. It is worth noticing the strong decrease ( Umean

Proper Orthogonal Decomposition
The flow field is decomposed for t = [2192s−2792s] within ±2R by applying POD to the velocity slices extracted 1R upstream for every 0.5s, a total of 1119 slices covering 10mins. Figure 2 shows the normalised POD eigenvalues and the cummulated eigenvalues, i.e. the amount of turbulent kinetic energy captured by the POD modes. The first four POD modes contains 22.6%, 6.8%, 5.6%, and 4.7% of the turbulent kinetic energy, approximately 40% of the total turbulent kinetic energy. The first 10 POD modes holds 55% of the turbulent kinetic energy, while it is necessary to include 215 modes to retain 90% of the total turbulent kinetic energy, corresponding to 0.1% and 19% of all POD modes respectively. The relatively large number of POD modes and mild slope is a clear indication that although the flow contains some larger structures, a large part of the turbulent kinetic energy stems from high frequency and small scale turbulence, which is distributed across a large number of POD modes.  Figure 3 shows the spatial POD modes 1,2,3, and 4 for all three velocity components. The first modes of U show clear dipole structures with extrema located outside the rotor plane. The first mode for V and W has maximum within the rotor area. The spatial modes for V and W are otherwise less distinct than those of U , and particular the second mode for V shows small and more scattered structures.
The temporal development or eigenfunction for POD mode 1 is given in Figure 4 along with the spectrum. A window was applied to the temporal eigenfunction before constructing the spectrum. Clearly, there are dominating frequencies, and the aim in the following section is to optimise the reconstruction by identifying and scaling these dominanting frequencies. Table 1 gives an overview of the optimised Reduced Order Model for the top four POD modes. The dominant frequencies for each POD mode is given along with the linear least squares optimised coefficients which gives rise to a coefficient of determination, i.e. the goodness of fit between the flow field reconstructed using POD( u) and with the Reduced Order Model( u). The goodness of fit is measured using coefficient of determination defined as

Reduced Order Model
where u 0 is the mean of u. Finally, the fraction of turbulent kinetic energy captured by the POD reconstruction is given.   Using only five distinct frequencies, it is possible to reconstruct 22% of the total turbulent kinetic energy with R 2 = 77%, where the scaled top frequency f = 0.00832Hz accounts for 48% of the total energy in POD mode 1. The coefficient of determination naturally drops each time another POD mode is introduced, since more energy is introduced into the reconstruction. The influence of the otherwise dominant frequencies on the flow for higher POD modes essentially disappears for POD modes higher than 4, i.e. it is not possible to truncate the higher modes in terms of a few dominant frequencies. This clearly demonstrates how the turbulent kinetic energy is spread over a large number of POD modes and that the large turbulent structures are captured by the first four POD modes. Therefore, only the top 4 POD modes and corresponding frequencies are used in the following to reconstruct the flow. The table gives the top nine frequencies, which reconstructs 40% of the turbulent kinetic energy contained by these large scale structures with R 2 = 0.60. Reconstructing using the dominant 18 frequencies would give a goodness of fit of 69%, but it is a tradeoff to balance the aim of reconstructing the overall large scale dynamics with the fewest possible number of variables. Therefore, only the 9 most dominant frequencies and top 4 POD modes are used in the following reconstruction, which corresponds to using only 0.4% of the frequencies for the full POD reconstruction.
Furthermore, it is noteworthy that several of the dominant frequencies for the individual   This indicates that the large scale dynamics of the flow originates from combining these spatial dynamics using similar frequencies.  Figure 5 compares the streamwise velocity from CFD(black), POD(red) and the Reduced Order Model(blue) reconstruction at the rotor centerline. As seen, there are fast and large fluctuations which are not captured by neither the POD nor the Reduced Order Model reconstructions. However, both clearly capture the low frequent dynamics. This is also shown in Figure 6, which shows snapshots in time of the streamwise velocity as given by CFD and the Reduced Order Model. The wake is very dynamic and although the Reduced Order Model only includes a limited number of frequencies and POD modes, it generally captures both the location of the wake deficit and the distortion of the wake, hence the large scale dynamics of

Discussion and Conclusion
The wake interaction and performance of a wind turbine in the interior of a wind farm is investigated using numerical simulations of an idealised case of an infinitely long row of wind turbines by applying cyclic boundary conditions. The wind turbines are modelled using the Actuator Line method coupled with the aeroelastic code Flex5, which calculates loads and performance of a NM80 wind turbine. Large Eddy Simulation is employed to model the subgrid scale turbulence. A Reduced Order Model is constructed based on a Proper Orthogonal Decomposition of the flow. POD yields the large scale spatial structures as well as temporal eigenfunctions. The Reduced Order Model consist of combining the spatial structures with a truncated representation of the temporal eigenfunctions by choosing the dominant frequencies for each POD mode. The Reduced Order Model is shown to reconstruct the flow in a physically correct manner by capturing the low frequency dynamics using only nine frequencies distributed on the first four POD modes, which constitutes 40% of the turbulent kinetic energy. The reconstruction of the flow could potentially be combined with high frequency turbulence to account for the energy deficit and the more extreme fluctuations.