Investigation of the effect of inflow turbulence on vertical axis wind turbine wakes

The aerodynamics of Vertical Axis Wind Turbines (VAWTs) is inherently unsteady, which leads to vorticity shedding mechanisms due to both the lift distribution along the blade and its time evolution. In this paper, we perform large-scale, fine-resolution Large Eddy Simulations of the flow past Vertical Axis Wind Turbines by means of a state-of-the-art Vortex Particle-Mesh (VPM) method combined with immersed lifting lines. Inflow turbulence with a prescribed turbulence intensity (TI) is injected at the inlet of the simulation either from a precomputed synthetic turbulence field obtained using the Mann algorithm [1] or generated on the-fly using time-correlated synthetic velocity planes. The wake of a standard, medium-solidity, H-shaped machine is simulated for several TI levels. The complex wake development is captured in details and over long distances: from the blades to the near wake coherent vortices, then through the transitional ones to the fully developed turbulent far wake. Mean flow and turbulence statistics are computed over more than 10 diameters downstream of the machine. The sensitivity of the wake topology and decay to the TI and to the operating conditions is then assessed.


Introduction
The operation of Vertical Axis Wind Turbines (VAWTs) essentially involves unsteady aerodynamics, and thus the shedding of complex vorticity structures in their wake. Indeed, VAWTs exhibit a wake topology that is governed by both the spanwise lift variations (like their horizontal axis counterparts (HAWTs)) and the temporal lift variations that arise over the blade rotation. This leads to a strongly three-dimensional vorticity field in the wake of VAWTs, an aspect that could translate into a more intense wake decay. Another interesting aspect of VAWTs is that they are naturally indifferent to variations in the wind direction; this could clearly translate into a lower sensitivity to ambient turbulence. These characteristics of VAWTs and their wakes hint at a possible lead over HAWTs once they will be installed in a farm; this has, in particular, led to predictions of higher land-specific power densities [2,3].
Such promises, together with potential operational gains (maintenance costs, the disappearance of yawing actuation), have spurred some definite research momentum in VAWT aerodynamics, in the shape of experimental [4,5] and numerical [6,7] studies. However, because of their unsteady aerodynamics, VAWT simulation and modeling tools have not reached yet the maturity of those for HAWTs, e.g. the Blade Element Momentum (BEM) method. The wake phenomenon for VAWTs has only been approached recently, either from an experimental [5] or a numerical perspective [8,9]. Among those, very few works have been dedicated to the investigation of the large scale development of the wake [9], and none have thus far considered the operation of the machine within a turbulent wind. This paper addresses this very last point and follows directly the work of [9]. That work concerned a standard, medium-solidity, H-shaped machine; all simulations were performed without ambient turbulence, and rather considered the details of the vortex dynamics at work in the near wake and their influence on the wake decay. In the present work, we introduce a turbulent wind and perform a parametric study of its influence. In particular, this study highlights the behavior of the main flow structures identified in the absence of turbulence [9] and the relative contributions of the intrinsic vortex dynamics and the interplay between the wake vorticity and the ambient turbulence. We consider the same baseline VAWT geometry and keep it operating at its optimum Tip Speed Ratio (TSR); we vary the Turbulence Intensity (TI) and the nature of the turbulence: isotropic or anisotropic. The present work relies on the same numerical tools as in [9]: the large-scale and fine-resolution Large Eddy Simulations are performed by a state-of-the-art Vortex Particle-Mesh (VPM) method combined with immersed lifting lines [10]. We add the management of a turbulent inflow performed through the injection of a vorticity field in a Lagrangian fashion; this field is either pre-generated through the Mann algorithm [1] or is generated on-the-fly during the simulation.
This paper is structured as follows. We only briefly recall the Vortex Particle-Mesh (VPM) method with immersed lifting lines and turbulent inflows in Section 2. In Section 3, we study the effects of turbulence on the reference VAWT and also perform a verification of the synthetic turbulence approach. We close this paper in Section 4 with our conclusions and perspectives.
where ν is the kinematic viscosity, and T M is the Sub-Grid Scale (SGS) model, here a multiscale model solely acting on the small scales of the LES field [11]. The advection term of this equation is treated with particles while their right-hand sides are evaluated efficiently on a mesh. The Poisson solver for velocities operates in Fourier space and it simultaneously allows for unbounded directions and inlet/outlet boundaries [12]. The resulting hybrid method requires back-and-forth interpolations between the particles and the grid but does not detract from the good numerical accuracy (in terms of diffusion and dispersion errors) and the stability properties of a particle method. We refer to [13,14,11] for detailed descriptions. The generation of vorticity along the blades is accounted for through an immersed lifting line approach [10]. The approach is very much akin to a Vortex Lattice method but only over the span of a time step. It relies on the Kutta-Joukowski theorem (see e.g. [15]) that relates the developed lift L to the relative flow V rel and the circulation bound around the local 2D airfoil L = ρV rel × Γ and it discretizes the bound circulation and the shed vortex sheets with particles, which are then added to the bulk flow. In this work, we account for unsteady effects through a Leishman-Beddoes dynamic stall model [16], following the indications of Dyachuk [17] and Scheurich [6] who present coefficients for various airfoils validated in the particular case of VAWT.

Turbulent inflow
In order to model a realistic incoming inflow, we generate synthetic turbulence at the domain inlet. A perturbation vorticity field ω is injected at the inflow as particles. This Lagrangian treatment, as opposed to a boundary condition, is consistent with the rest of the methodology and does not impose any additional time step constraint: a several grid points-thick slab of turbulent wind can enter the domain over a single time step. This technique, originally developed in [10,18], is similar to that applied in other recent Vortex Particle methods works [19,20,21] but has been further developed here for more flexibility.
In a refinement of what was done in [18], the introduction of ω into the computational domain now relies on a smooth blending of ω with the domain vorticity, instead of a direct insertion of the ω particles into the domain. More specifically, a buffer computational domain co-exists with the main computational domain; this turbulent inflow contains particles which are progressively let free: their velocities and their stretching terms transition smoothly from a uniform flow at the prescribed inflow velocity U ∞ into the fully three-dimensional flow at the end of the buffer region. This change allows to handle arbitrary remeshing frequencies and non-uniform inflows.
In this paper, two ω generation techniques are used and compared. A first one relies on simply computing ω = ∇×u of a pre-generated database of velocities u , here using Mann's [22] algorithm. We introduce a second technique that avoids the pre-generation and storage of this database. At each time step, we generate a two dimensional slab of a turbulent velocity field using Mann's algorithm again; this velocity slab is correlated in time with the previous slabs using an asymmetric time filter as in [23] in order to produce realistic correlations in the streamwise direction.
Both methods produce fields that are periodic in the transverse directions; those are made compact through a smooth clipping, as we are working with unbounded transverse directions. Preliminary results for the second technique are provided below and are compared to those obtained using the pre-computed turbulence.

Parameters
This whole section focuses on the same low solidity H-VAWT studied numerically in [24,6], and then in [9] in configurations without any turbulence. For the sake of completeness, we here briefly recall its main parameters: an aspect ratio AR = H/D = 1.5, a solidity σ = nc/D = 0.1725, and constant-chord NACA0015 airfoils. We here only consider its optimum power operating point, TSR = ΩR/U ∞ = 3.21, as found by [25]. Unless otherwise specified, the numerical simulations involved D/h = 96 mesh points/particles per rotor diameter and extended up to 17 diameters downstream of the rotor axis. This discretization and the chosen turbulent inflow sizes lead to simulations with 3 10 8 to 7 10 8 particles/mesh points; those were run on 2560 CPU cores for wall-clock times between 30 (low TI cases) and 54 (high TI anisotropic case) hours. Energy (TKE), and of the anisotropy of this inflow, as characterized by the parameter Γ which is equal to zero for isotropic turbulence. The runs with isotropic turbulence I2 and I7.5 compare with wind tunnel conditions at low and intermediate grid-generated turbulence intensities; run A7.5, is a scenario with anisotropic turbulence at the same TI as I7.5, which is representative of a small-to-medium VAWT (D 20 − 50m) in atmospheric turbulence. All these turbulent inflows were pre-generated using Mann's algorithm [1].

Aerodynamics
We first consider the effect of these turbulence scenarios on the Angle-of-Attack (AoA) and on the normal and tangential force coefficients at mid-height, respectively F n /(c q 0 ) and F t /(c q 0 ), where q 0 is the dynamic pressure of the mean free stream. We will take, as reference, the behavior without ambient turbulence computed in [9]. We recall here the main features of these quantities over a blade rotation in the absence of turbulence: (1) the behavior is not symmetrical for the upstream and downstream legs because of the reduced velocity encountered downstream; (2) the downstream leg also exhibits oscillations in the AoA and force coefficients as the blade encounters vortex sheets shed during the upstream leg. (c) Tangential force coefficient Figure 1. H-type VAWT aerodynamics: profiles of the angle of attack and of the normal, F n /(c q 0 ), and tangential, F t /(c q 0 ), force coefficients at mid-height versus the blade angular position θ for the turbulent inflows: I2 (----), I7.5 (-· -), and A7.5 (· · · · · ·) and the reference without turbulence [9] (--); the two-standard deviations (±2σ) envelopes are indicated by the colored areas.
The average aerodynamic behaviors do not deviate markedly from the reference case without turbulence. The downstream oscillations decrease with the TI, as the blades interact with the upstream-shed structures at more randomized locations. The standard deviations globally exhibit orders of magnitude consistent with the turbulence level. They are however modulated by the blade position: the angle-of-attack is perturbed the most when the blade travels perpendicularly to the freestream. The angles-of-attack are practically zero around θ = 20 • and 180 • ; the zero slope of the drag curve at α = 0 then translates such values into a collapse of the tangential forces perturbations, i.e. very small envelopes (Fig. 1(c)). Finally, we note the strong difference between the deviations achieved in the isotropic and anisotropic cases. For equal TIs, an anisotropic inflow redistributes more of the turbulent kinetic energy to the velocity components which impact the aerodynamics the most: those in the horizontal plane (u 1 and u 3 ).
The mean power coefficient C p = P/( 1 2 ρAU 3 ∞ ) obtained without inflow turbulence was C p = 0.339 [9]. When inflow turbulence is added, the obtained power coefficients are 0.346, 0.342 and 0.357 for the cases I2, I7.5 and A7. 5 without turbulence is mainly due to the increased averaged energy flux in the turbulent cases: these cases maintain the averaged inflow velocity atū = U ∞ but the average u(u 2 + v 2 + w 2 ) differs from U i nf ty 3 . Yet, the C p of case I7.5 is slightly lower than that of case I2, denoting a small decrease of the aerodynamic efficiency due to the larger fluctuations.
3.3. Wakes 3.3.1. Vortex dynamics Instantaneous visualizations of the wakes are shown in Fig. 2 for the considered TIs. The turbulence-free baseline (Fig. 2(a)) provides a clear view of the fundamental vortex dynamics at work in the near wake. On the sides, the blade sheds vortex sheets due to the temporal variation of the circulation, dΓ/dt, and at the top and bottom, it sheds tip vortices with varying circulations. The near-wake is then governed by the instabilities and the merging phenomena of vortices of unequal circulations, all of which can be seen to originate in the corners of the wake. We refer to [9] for a more thorough discussion. The ambient turbulence affects the wake at several scales. In the near-wake, one can observe the disturbances on the shed tip vortices (see the zooms in Fig. 2(b) to 2(d)) increasing with the TI. From very slight, they increase markedly for the high TI cases, and even lead to reconnections between successive tip vortices in the vicinity of the blades in the downwind part of their rotation. From a global point of view, if an absence of turbulence produces a quite periodic array of vortices and reconnection events in the corners of the wake, increasing the turbulence level leads to more irregularly spaced and distorted vortices and reconnections occurring on the top and bottom of the near-wake. One thus expects a more distributed production of turbulence in the near wake.
The wake behavior at large scales exhibits meandering with amplitudes increasing with the TI, and we note that it is already present in the baseline without turbulence. The wavelength of this meandering motion is of the order of 5D and thus much larger than the integral length scale of the synthetic turbulence boxes (D/3 and D). We leave the detailed investigation of this meandering behavior as a subject of future work.

Statistics
We study the average behavior of the wakes by means of the mean axial velocityū, and the mean turbulent kinetic energyk; these statistics were collected over a period T avg = 24 D/U ∞ , corresponding to the flow-through time of the pre-computed inflow turbulence field. Figure 3 shows a horizontal slice of these statistics for the studied TIs and for the baseline. The previous study had identified a backflow region in the wake (Fig. 4(a)) at intermediate and high TSRs. This feature is here shown to persist for the isotropic turbulent inflows; it also moves closer to the rotor as the TI is increased, going from 6 D for no TI to 4 D and less for TI = 2% and 7.5%. This structure disappears for the anisotropic TI case; its large scale perturbations are then able to disrupt the near-wake. The TKE signatures translate and quantify the above-mentioned accelerated onset of vortical instabilities in the near-wake. Even at a low TI, TKE production increases on the sides of the wake and progresses upstream towards to the rotor (Figs. 4(b) and  3(c)). The TKE of the anisotropic case exhibits a strongly increased intensity; its proximity to the rotor and the angular spread of its signature are indicative of more pronounced large scale wind direction changes.
This higher dissipation translates in an increased smearing of the velocity deficit, with the most important gain for the anisotropic case; this is confirmed in the transverse slices of Fig. 4. We note the more even distribution of TKE in the periphery of the near wakes with a slight early bias for the blade-travelling-downwind side. The bean-shaped distributions for the deficit and decay are well preserved for the isotropic cases, up to 10 diameters. For the anisotropic case, the wake has already much decayed at 10D downstream.

Preliminary results on turbulent inflow comparison
The generation of a turbulent inflow using a precursor has two drawbacks. First, it requires to generate a 3D turbulence database, store it and manage its input into the code. Second, as this database has a finite streamwise extension, the information entering the flow domain is repeated in time leading to periodic repetition of flow structures. With the alternative method  investigated here, these two drawbacks are overcome as information is generated on the fly. Comparisons of flow diagnostics (T I = 7.5%) obtained using the two methods are provided in Figs 5 and 6: a good agreement is found between both approaches. This comparison constitutes a first step towards validation of the new proposed approach.

Conclusions
A Vortex Particle-Mesh method with immersed lifting lines, here briefly presented, has been applied to large scale and high resolution LES of VAWT wakes with various turbulent inflows. The method is capable of tracking vortical structures over very long times and distances. Two TI levels (2.5% and 7.5%) have been considered for isotropic inflow turbulence whereas anisotropic turbulence has also been considered for the higher TI case. For each case, the mean flow and mean TKE fields have been obtained. Increasing the turbulence level accelerates the transition of the wake and this tends to bring the low speed recirculating region closer to the turbine. For the higher TI considered (7.5%), the anisotropy of the turbulence has a strong effect since the decay is much faster in the anisotropic case than in the isotropic one. The recirculation zone is also suppressed in the case of the anisotropic turbulence and the TKE in the wake is also much higher. The wake meandering is also found to increase with the TI level, and to be further enhanced by the turbulence anisotropy. These large numerical databases contain rich information that will be further analyzed.