Turbulent large-scale structure effects on wake meandering

This work studies effects of large-scale turbulent structures on wake meandering using Large Eddy Simulations (LES) over an actuator disk. Other potential source of wake meandering such as the instablility mechanisms associated with tip vortices are not treated in this study. A crucial element of the efficient, pragmatic and successful simulations of large-scale turbulent structures in Atmospheric Boundary Layer (ABL) is the generation of the stochastic turbulent atmospheric flow. This is an essential capability since one source of wake meandering is these large - larger than the turbine diameter - turbulent structures. The unsteady wind turbine wake in ABL is simulated using a combination of LES and actuator disk approaches. In order to dedicate the large majority of the available computing power in the wake, the ABL ground region of the flow is not part of the computational domain. Instead, mixed Dirichlet/Neumann boundary conditions are applied at all the computational surfaces except at the outlet. Prescribed values for Dirichlet contribution of these boundary conditions are provided by a stochastic turbulent wind generator. This allows to simulate large-scale turbulent structures - larger than the computational domain - leading to an efficient simulation technique of wake meandering. Since the stochastic wind generator includes shear, the turbulence production is included in the analysis without the necessity of resolving the flow near the ground. The classical Smagorinsky sub-grid model is used. The resulting numerical methodology has been implemented in OpenFOAM. Comparisons with experimental measurements in porous-disk wakes have been undertaken, and the agreements are good. While temporal resolution in experimental measurements is high, the spatial resolution is often too low. LES numerical results provide a more complete spatial description of the flow. They tend to demonstrate that inflow low frequency content - or large- scale turbulent structures - is an important parameter when simulating wake meandering and plays a significant role.


Introduction
Wind turbine wakes are critical for the design and operation of wind farms [1]. The wake velocity deficit contributes to decrease the wind farm energy production [2,3]. Wakes should be considered during the design of wind farms in order to decrease energy losses associated with velocity deficit [4]. Wakes also increase turbulence intensities which influence the structural behavior and the life time of a wind turbine. Siting studies should consequently consider wakes during their analyses. Finally, the recovery rate of wake velocity deficit is highly dependent on the level of turbulence intensity [5]. Low turbulence intensities yield low recovery rate. This explains why wakes are critical for offshore developments.
Wakes can be considered as steady or unsteady phenomena. Analysing wake behavior as a steady phenomenon [6,7,8] typically leads to an overestimation of the velocity deficit. A more complete analysis is to consider wakes as non stationary flow structures [9,10,11,12]. One of the key aspect of wind turbine wake dynamics is the phenomenon of meandering [13,14,15] which contibutes to fatigue loading. Measurements and observations of this phenomenon have been undertaken in situ [13,15] and in wind tunnel [16,17]. The causes of wake meandering were associated experimentally with the unsteady behavior of large-scale turbulent structures [18,19] and the tip vortex instabilities in turbulent wake flows [20]. A pragmatic and efficient wake meandering model based on the effects of large scale (low frequency) structures has been proposed by Larsen et al. [10]. España et al. [18] have shown experimentally the importance of large scale turbulent structures on wake meandering behind porous disks.
When developing a simulation technique for wake meandering, an important aspect to consider is the flow model, including turbulence modeling, and the way wind turbine action on the flow is modeled with its associated level of details. Regarding flow and turbulence models, it ranges from basic momentum conservation to direct simulation of Navier-Stokes equations. Regarding the turbines, an actuator disk or actuator lines can be considered.Various combinations have been proposed and successfully applied [21,22,23]. Each combination was proposed in order to efficiently simulate a specific phenomenon (e.g. tip vortex shedding, blade tip losses, wakes). References [23,24] are relevant to the work presented in this paper since they describe a simulation technique specifically developed to simulate unsteady wake behavior.
In this work, a simulation technique suitable for the analysis of large scale turbulent structures in Atmospheric Boundary Layer (ABL) is proposed and presented. A crucial element of this simulation technique is the generation of the stochastic turbulent atmospheric flow. This is an essential capability since one supposed source of wake meandering is these large -larger than the turbine diameter -turbulent structures.The unsteady wind turbine wake in ABL is simulated using Large Eddy Simulations (LES) in which the wind turbine is represented by an actuator disk. In order to dedicate the large majority of the available computing power in the wake, the ABL ground region of the flow is not part of the computational domain. Instead, mixed Dirichlet/Neumann boundary conditions are applied at all the computational boundary surfaces except at the outlet. Prescribed values for Dirichlet contribution of these boundary conditions are provided by the stochastic turbulent wind generator. This allows to simulate large-scale turbulent structures -larger than the computational domain -leading to an efficient simulation technique of wake meandering. Since the stochastic wind generator includes shear, the turbulence production is included in the analysis without the necessity of resolving the flow near the ground. The classical Smagorinsky sub-grid model is used. It has been determined that in such an implementation, this choice of sub-grid model is not critical since more than 80% of the kinetic energy of turbulence is resolved.

Mathematical model 2.1. Turbine representation
The turbine is modelled using an actuator disk of diameter D with an axial induction factor a = 0.19 corresponding to the measured value presented in [17]. This disk extracts momentum from the flow which produces a local deceleration. This constitutes an acceptable approximation of the wind turbine far wake, as described by Froude's theory [25]. This disk is introduced in the form of an internal boundary condition at the location where it is inserted, and where a pressure discontinuity is prescribed according to Eq.(1) but through which the velocity is continuous.
The imposed pressure discontinuity ∆p is computed using the fluid density ρ, the instantaneous local wind speed at the disk location, U d and the prescribed induction factor a:

Governing equations and turbulence modelling
The governing equations are the Navier-Stokes equations in which gravity and Coriolis terms are neglected. The LES approach with implicit filtering is adopted in this work. Consequently, only large-scale turbulent structures are resolved, and the smaller scale structures are modelled using the classical Smagorinsky model with a model constant C s set to 0.167.

Synthetic turbulence generation
The stochastic generation of turbulent velocity series is a subject that has been approached by various authors [26,27,28], but the most commonly employed in the field of wind turbine simulation is probably that of Mann [29], which relies on the inversion of the spectral tensor of atmospheric turbulence. In this work, the turbulent inflow is obtained by mean of a stochastic method similar to that of the aformentioned authors, which is described in [30]. The resulting series present realistic turbulence characteristics, especially regarding the largest scales of turbulence, which makes it particularly useful for the simulation of atmospheric flows. It makes use of an inverse wavelet transform to enable the simulation of sheared, non-homogeneous turbulent fields. It also introduces the concept of motion compensation in order to incorporate the temporal dynamic of the production and decay of the turbulent eddies.

Numerical method
All simulations are carried out using OpenFOAM [31]. The gradient terms are discretized using second order linear interpolations. The divergence terms are also discretized using secondorder linear interpolations with the possibility of adding up to 20% of first-order upwinding to prevent any oscillations. This is the option filteredLinear in OpenFOAM. Time marching is done using backward Euler scheme. The pressure implicit split operator (PISO) algorithm [32] is used to handle pressure-velocity coupling. In OpenFOAM, it is available through pisoFoam. The systems of velocity equations are solved using linear algorithm (smoothSolver ) while the geometric-algebraic multigrid algorithm (GAMG) is used for pressure equations.
3.1. Boundary conditions Section 2.3 briefly described the methodology developed in order to produce realistic time series of turbulent velocity in ABL using a stochastic generator. In this section, the procedure to introduce these time series into the LES through appropriate boundary conditions is described.
There are two types of boundary conditions (BC) for differential equations: Dirichlet and Neumann types. For Dirichlet type of BC, the value of the variable is fixed at the boundary. For Neumann type of BC, it is the gradient of the variable that is fixed at the boundary.
Dirichlet and Neumann types of boundary conditions (BC) can be combined in a so-called hybrid BC. For a general dependent variable µ -which can be the pressure or velocity components -using f as the fraction of the BC associated with the Dirichlet's type of BC and the indices c and f to indicate the variable location at cell center and boundary face, respectively, the hydrid BC can be written as: Dirichlet BC, prescribed value of µ: µ Dirichlet = µ.
Neumann BC, prescribed ∂µ f /∂x , x being the position: Hybrid BC combining the previous two BCs: The approach adopted in [23,33,24] consists in using steady-state BCs and turbulence is introduced through unsteady source terms in a vertical plane lying upstream of the rotor in the computational domain.
An alternative approach to introduce the turbulence is adopted in this work. It consists in imposing the fluctuating velocities at all boundaries of the computational domain except at the outlet. The velocity values are extracted from a spatially and temporally coherent synthetic velocity field, which is generated according to the stochastic method described in Section 2.3. This approach ensures a reasonable unsteady behavior of the flow in the entire computational domain: inspection of the numerical solutions reveals that the resolved turbulent kinetic energy near the boundary is similar to that which can be observed further inside the domain. Additionaly, it makes it possible to introduce, and study, turbulent structures larger than the computational domain. The largest turbulent scales in the synthetic velocity field may be larger than the size of the LES computational domain. Since the largest turbulent scales of the synthetic field are applied at the inlet, lateral, top and bottom boundaries, it ensures that the computational domain is swept by these large turbulent structures, and consequently, the effect of these large structures on the meandering of the wake is reproduced in the computational domain. However, coherence between the turbulent structures in the computational domain and the external atmospheric turbulence is needed. This is verified by comparing the time-spectra of the synthetic velocity and resolved velocity inside the LES as shown on Fig.4. It shows that the turbulence levels are similar for most frequencies, and especially low frequencies (large scales). Since the turbulence levels are similar, it is reasonable to assume that the simulated velocity in the domain is directly driven by the synthetic field, with a slight amount of discrepancy depending on the scale of the structures. This discrepancy is caused by the fact that the evolution of the simulated velocity field is obtained by solving the Navier-Stokes equations, whereas the evolution of the synthetic velocity field is the result of a simplified advection scheme relying mainly on Taylor hypothesis.
Consequently, the synthetic velocities are prescribed on all boundary faces of the computational domain, except at the outlet face which is orthogonal to the main wind speed direction. Zero gradient is applied on velocities at the outlet. For the lateral, top and bottom faces, the average mass flow rate is zero on each boundary face, but since the flow is turbulent, instantaneous flow rate may be non-zero.
During the initial phase of the development, a Dirichlet type boundary condition (fixedValue in OpenFOAM) was imposed both at inlet and lateral faces. The prescribed velocity values were given by the stochastic generator. This implementation allowed to obtain overall physically meaningful solutions, but with significant spurious spatial oscillations near cells adjacent to boundary faces with instantaneous outflow. This undesirable behavior of the solution is produced mainly by an inconsistency between the small-scale structures computed by the LES and the stochastic generator: the discrepancy between the synthetic velocity field and the simulated LES velocity field is mainly produced by the advection which is calculated differently in LES and synthetic turbulence generator. Only a minor part of this discrepancy can be attrubuted to the difference in resolution between synthetic and simulated turbulence.
While large-scale structures are typically coherent between the solutions of LES and stochastic generator, it is not the case for small-scale structures which experience a different evolution since their coherent time-scales are much smaller.
This inconsistency can be solved by using hybrid boundary conditions where the Dirichlet fraction f in Eq.(2) is taking values from 1 -for pure Dirichlet -to 0 -for pure Neumanndepending on the sign of the local wall normal velocity component (and thus the flow through the boundary face) of the instantaneous flow as given by Eqs. (3) and (4). It is to be noted that there is a value of f for each velocity component. This hybrid boundary condition is applied to the components of the velocity tangential to the wall, and is based on the value of the wall- normal component of the velocity. Its role is to loosen the specified value of the velocity at the boundary during outflows in order to prevent spurious spatial oscillations. Regardless of the direction of the flow, the wall-normal component of the velocity is always specified with a pure Dirichlet condition. When the normal component was not specified, a decrease of theaverage velocity of the flow in the main direction was observed due to an excess of mass exiting the computational domain. Following this reasoning, Eq. (2) can be re-written in a vectorial form by replacing the scalar quantity f by the tensor F . Using tensor F allows the specification of a different value of the Dirichlet fraction for each velocity component expressed in term of the face normal n oriented toward the interior of the computational domain: The operator ∧ is the external product, resulting in a tensor. I is the identity matrix. Eq.(2) can be written as: The values of u Dirichlet at all boundary faces are given by the stochastic generator at proper time and location. Values of u N eumann at all boundary faces are determined from the velocity gradient computed from the averaged flow field as given by the ABL logarithmic expression. This hybrid boundary condition allows to maintain the flow rates at all boundaries as those given by the stochastic generator. Zero gradient is set at boundaries for pressure with a reference pressure prescribed at a given point in the computational domain.

Computational domain
The computational domain is illustrated in Fig. 1. The bottom face is not located at the ground, but at a height of 5m where turbulence is assumed to be completely developed. As such, the velocity is not zero at this bottom face. Coarser mesh resolutions (5.66m and 6.82m) were also tested. The results were very close. Furthermore, given the simplicity of the WT model used here, a finer mesh does not appear necessary. While we do focus on the large turbulent scales and do not make the wake recovery and diffusion the main priority of this study, it is important to point out that the values of the velocity on the centerline of the wake at 2D and 5D behind the rotor are in very good agreement with experimental data, as shown in Section 4.2.3.
Computational domain parameters such as domain width Ly, domain height Lz, rotor diameter D, number of grid points N grid , grid resolution ∆ are reported in Table 1. The relatively small domain height adopted is the result of a compromise between computational cost and modelling accuracy. The choice of a relatively small domain height has been made to allow the use a uniform mesh to ensure that the implicit turbulence filtering scale does not change across the domain. Given such a constraint for the mesh, increasing the domain size while keeping the same resolution quickly increases the number of cells. Moreover it was decided to use a domain width greater than height so as to prioritize horizontal wake meandering over vertical, since the horizontal meandering shows a greater amplitude.

Actuator disk implementation
The internal face at the disk location is determined by selecting a vertical plane in the grid in order to ensure that the disk is composed of existing grid cell faces. A grid face is considered to be part of the actuator disk when its center lies within a circle having the same diameter as the actuator disk. All grid faces having these caracteristics compose the discretized actuator disk on which the internal boundary condition is applied. Since the grid is cartesian, the discretized actuator disk is not a perfect circle. Its perimeter is of stair shape. With a grid resolution of 5m, there are about 16 cells radially over 1D. While this might seem coarse, an uniform Cartesian mesh is simple and has the advantage of ensuring that the implicit LES filter scale is homogeneous accross the whole domain. Through the actuator disk, only the pressure is discontinuous, the velocity is continuous. This does not create a singularity in the code, as it is implemented as a pair of coinciding cyclic boundaries internal to the domain. In order to allow comparisons with experimental measurements in wind tunnel [17], a disk of 80 m immersed in an ABL having a roughness height z 0 = 0.02m is analysed. The actuator disk is located downstream of the inlet computational domain at a distance of 4D. This distance is sufficient for the turbulence to reach an equilibrium state. Figs.2 and 3 shows the vertical profiles of velocity and turbulent intensity at various longitudinal locations downstream of the inlet. The horizontal homogeneity is almost perfectly obtained, indicating that the proposed top and bottom BCs behavior is adequate. There is a noticable, but small, acceleration near the ground and a deceleration near the top region of the computational domain. Transition from acceleration toward deceleration occurs at about one third of the computational domain height.

ABL simulations
The frequency weighted longitudinal velocity energy spectra in Fig.4 shows the presence of large-scale structures. The theoretical slope of −2/3 is also observed at high frequencies with high-frequency deviations most likely caused by contamination from the LES sub-scale filtering.    Their amplitude was of the order of 10% the average flow velocity, so not particularly large but definitely unphysical. Such spurious oscillations would clearly be seen on instantaneous as well as average velocity contours.

Wake simulations 4.2.1. Mean profiles
The vertical profiles of the longitudinal velocity, Fig. 6 are consistent with the ABL velocity profiles previously presented. They exhibit an increase in velocity in the profile outside of the WT wake, which is caused mainly by the relatively small height of the computational size with respect to the disk diameter. This choice of computational domain height is however deemed acceptable, since the main goal of this work is to study the wake meandering. This decision is further supported by the results of a study of the sensitivity to the domain size and to the mesh resolution. The results of this study are reported in [30].
The values of velocity at the first cell near the ground at position z = 5 m are almost identical for all streamwise locations. This is due to the implemented boundary conditions, which let the tangential velocity be computed from the LES solution when the flow is leaving the computational domain, but set the velocity vector to values prescribed by the stochastic generator when the flow is entering. The value of the longitudinal velocity is effectively fixed 50 % of the time. In fact the probability density of the synthetic velocity normal component at the   The implemented porous disk modeled is having a non-uniform pressure jump calculated from instantaneous local velocity. This is clearly not identical to an actuator disk having a uniform pressure. The current implementation is expected to be closer to an actual turbine, in that the average wind speed varies with height which produces a change in the average pressure jump. Fig. 7 presents the vertical profile of the streamwise turbulence intensity at various downstream locations. The actuator disk is located at x = 0D. The value of the turbulence intensity in the middle part of the near wake is very close to the upstream undisturbed ABL value. This is to be expected since the main turbulence production mechanism is associated with velocity gradients which are significant only at the wake edges in the near wake region.

Meandering
It is important to note that meandering is considered in this work as the time variation of the location of the wake velocity deficit centre caused by the large scale structures in the incident ABL. The unsteady simulation of the actuator disk immersed in the ABL allows the observation of wake time evolution. Fig. 9 presents the instantaneous flowfield in a vertical plane in which the turbine axis of rotation lies. It shows a variation of the wake trajectory which is typically associated with meandering as defined previously. The procedure adopted here to determine the instantaneous position of the wake is the same as the one described in Ref. [17]. It is different than the 2D Gaussian fitting procedure described in [19] but should yield similar results. It allows the calculation of wake position temporal series along transversal coordinates. The effets of averaged shear on velocities are removed when calculating z wake position in order to remove the effects of the undisturbed ABL velocity profile. Fig. 8 shows coherence between the transverse velocity measured at the center of the actuator disk and the instantaneous location of the wake centerline on the y axis for different distances d on the x axis downstream of the rotor. The usual methods to improve the statistics have  Figure 9. Wake instantaneous velocity field in median vertical plane for a = 0.19 and z 0 = 0.02 m been used where appropriate. The main reason why the spectra of Fig. 8 are noisy is that there are no symmetries in the flow that can be exploited to perform spatial averaging of the turbulence spectra measured in the wake of the turbine. Spectra from Fig. 8 looks noisier than those shown in Fig. 4. This is mainly associated the the scales which are of log type in Fig. 4 and linear in Fig. 8. The behavior of computed coherence shown in Fig.8 is very similar to the one computed from the observations [17] with high values of coherence at low frequencies and values close to zero at high frequencies. This change in behavior is occuring at a nondimensional frequency around of 0.3 − 0.4 (corresponding to a length scale of less than 3D), which is again in agreement with the wind-tunnel observations [17]. This frequency correponds to a wavelength of about three times the diameter of the rotor. The coherence is decreasing with the dowstream distance from the disk. The curves presented in Fig.8 are more noisy than the coherence obtained from the measurements. This is attributed to the reduced total number of computed samples/datapoints/timesteps in comparison with the experimental measurements. A better statistical convergence could certaintly be achieved using simulations with longer time periods, but considering the non-linear convergence rate with respect to the number of samples/datapoints/timesteps, the cost of these simulations becomes rapidly prohibitive. The physical time simulated to obtain the coherence is approximately 1 hr, which is already a long period in turbulent ABL simulations.

Validation
The numerical results have been compared to experimental observations presented in Ref. [17]. More specifically, the following characteristics are compared: • U 2D /U ∞ : ratio between the average wake centre velocity 2D dowstream and the upstream undisturbed velocity at hub height. • U 5D /U ∞ : ratio between the average wake centre velocity 5D dowstream and theupstream undisturbed velocity at hub height. • σ wake,y /D : standard deviation of the wake centreline horizontal location at 5D downstream.
• σ wake,z /D : standard deviation of the wake centreline vertical location at 5D downstream.
Numerical and experimental values are reported in Table 2.
Regarding non-dimensional velocities, the agreement with the measurements is almost perfect. However, the error on the standard deviation of the wake centre horizontal position is more than 20%. Nevertheless, these results shows that the meandering is qualitatively reproduced. The error on the vertical meandering is much larger. It is expected that some improvement  Table 2. Experimental/numerical comparisons could be brought to these simulations by improving the synthetic turbulence generator: the spectral tensor model used in this work quite severely underestimates the fluctuations of the vertical components of the velocity, which we think is responsible for the large discrepancy found for the standard deviation of the vertical meandering. However, the integral scales from the numerical simulations and the measurements are similar (around 200m at 2.5D). It is therefore suspected that there is a disagreement between the numerical and experimental turbulent structure lifespans. One of the source of discrepancy is suspected to be the characteristics of the incident ABL turbulence. While the turbulence intensities are identical between experiments and simulations, the spectral tensor implemented in the stochastic generator is relatively simple and consequently does not reproduce some anisotropic characteristics of the wind tunnel turbulence. Using the spectral tensor of Mann [29] may improve the results of the simulations. Despite these limitations, the results shown are still useful given that there are relatively few actual results available in the literature on the subject of wake meandering. Moreover, it is believed that the proposed methodology is physically sound and that the current limitations can all be addressable in the near future.

Conclusion
This paper presents an application of stochastic temporal velocity series for the LES simulation of wind turbine wakes in atmospheric boundary layer. The actuator disk is modeled by means of a pressure jump dependent on the local unsteady velocity and the prescribed induction factor, which is the similar numerical version of the porous disk used experimentally.
The fundamental role of boundary conditions based on a directional hybrid formulation, which adopts the behavior of a Neumann or Dirichlet type boundary conditions depending on the direction of the flow at the boundaries, has been demonstrated. It is shown that this type of BC ensures at the same time a strict control of the mass flow rate in the computational domain, while preventing spurious spatial oscillations near the boundaries with locally exiting flow.
The study results show that the method is effective in producing longitudinally homogeneous velocity profiles in the entire computational domain. The turbulent intensity profiles are also exhibiting a relatively fully developed behavior.
The proposed method allows to model realistic atmospheric flows using LES, without requiring a precursor simulation which can be costly in term of computing resources and challenging to implement. Such a method could be useful in many engineering applications in the ABL.
The instantaneous horizontal or vertical position of the wake on a transverse section is determined from the instantaneous velocity deficit profile whose values are used to calculate a center of gravity of the position of the wake. The result is a time series of the transverse position of the wake center. Coherence between the wake velocity time series and the upstream transverse velocity reveals a very strong correlation of low frequency which decreases with increasing frequency and reaches a non-significant value for any wavelength of less than about two times the diameter of the rotor.
The relatively large values of the measured coherence suggest that the real-time forecasting of meandering from upstream transverse velocities is a realistic prospect which could be done, for example, from global stress measurements on the mast or wind turbine rotor.
The comparison with experimental results shows that the computed velocity deficits are quite satisfactory. The simulated horizontal meandering qualitatively reproduces the observed experimental behavior. The agreement of the numerical vertical meandering is less convincing. Improved synthetic turbulence from the stochastic generator should increase the agreement.