Large Eddy Simulation of Vertical Axis Wind Turbine Wakes

In this study, large eddy simulation (LES) is combined with a turbine model to investigate the wake behind a vertical-axis wind turbine (VAWT) in a three-dimensional turbulent flow. Two methods are used to model the subgrid-scale (SGS) stresses: (a) the Smagorinsky model; and (b) the modulated gradient model. To parameterize the effects of the VAWT on the flow, two VAWT models are developed: (a) the actuator swept-surface model (ASSM), in which the time-averaged turbine-induced forces are distributed on a surface swept by the turbine blades, i.e., the actuator swept surface; and (b) the actuator line model (ALM), in which the instantaneous blade forces are only spatially distributed on lines representing the blades, i.e., the actuator lines. This is the first time that LES has been applied and validated for the simulation of VAWT wakes by using either the ASSM or the ALM techniques. In both models, blade-element theory is used to calculate the lift and drag forces on the blades. The results are compared with flow measurements in the wake of a model straight-bladed VAWT, carried out in the Institute de Méchanique et Statistique de la Turbulence (IMST) water channel. Different combinations of SGS models with VAWT models are studied, and a fairly good overall agreement between simulation results and measurement data is observed. In general, the ALM is found to better capture the unsteady-periodic nature of the wake and shows a better agreement with the experimental data compared with the ASSM. The modulated gradient model is also found to be a more reliable SGS stress modeling technique, compared with the Smagorinsky model, and it yields reasonable predictions of the mean flow and turbulence characteristics of a VAWT wake using its theoretically-determined model coefficient.


Introduction
Revived attention towards vertical axis wind turbines (VAWTs) has been observed in recent years.VAWTs are claimed to have several advantages over the conventional horizontal axis wind turbines (HAWTs), such as: insensitivity to the wind direction and, consequently, the absence of any yaw equipment; lower noise, due to the lower tip-speed of the blades; and placement of the drive train system on the ground [1].Recently, Dabiri [2] has claimed that the power density of wind farms consisting of counter-rotating VAWTs can be potentially greater than HAWT wind farms by an order of magnitude.To understand the performance of VAWTs better, detailed knowledge of the flow through them and the wakes behind them is necessary.Furthermore, to optimize the efficiency of a single VAWT or an array of VAWTs in a wind farm, an accurate and reliable prediction of the flow field behind the turbines is essential.For this purpose, numerical simulations can provide valuable quantitative understanding of VAWT wakes in a relatively fast and cheap way.
Almost 40 years after Darrieus' patent was registered (in 1931), researchers, mainly from some North American institutes (such as the National Research Council of Canada, NASA Langley Research Center and the Sandia National Laboratories), started to study VAWTs extensively.Rigorous and exhaustive experimental data from both wind tunnel and field measurements for different sizes of VAWTs were collected by these institutes in the 1970s and 1980s [3][4][5][6][7][8].These experiments and other wind tunnel and field measurements, which were done during those years [9][10][11][12][13][14], mainly focused on overall rotor performance (e.g., power and torque) and loading on the blades.The first attempts to experimentally study the wake behind a VAWT date back to the late 1970s, when Muraca and Guillotte [15] performed velocity measurements in the near wake of a VAWT in a wind tunnel, and Vermeulen et al. [16] carried out flow measurements in the wake of a full-scale VAWT in the natural wind environment.None of them reported wake turbulence data.Later, in the 1980s, Strickland et al. [17] and Brochier et al. [18] carried out water channel experiments to quantify the flow field downstream of a model straight-bladed VAWT.Bergeles et al. [19] performed hot-wire measurements in the near wake of a Troposkein-shaped VAWT in a wind tunnel and reported the profiles of mean velocity and turbulence intensity for different tipspeed ratios and two downstream positions.More recently, Battisti et al. [20] carried out wind tunnel measurements on the wake of a VAWT and provided data on time-averaged and phase-resolved velocity and turbulence intensity in a vertical plane located 1.5 rotor diameters downstream of the turbine axis.They showed that large-scale aerodynamic unsteadiness is greater in the tip region than in the central region.With the emergence of particle image velocimetry (PIV) techniques, some flow visualizations and measurements have been done to study the dynamic stall phenomenon and vortex evolution inside and behind a VAWT [21][22][23][24].Tescione et al. [24] has reported contours of average velocity and vorticity up to about 1.5 rotor diameters downstream of the turbine in a wind tunnel.
Besides experimental investigations, analytical and numerical studies have also been employed with the aim of VAWT performance prediction.The aerodynamic models that have been proposed for VAWTs can be classified in two main groups: streamtube models [1] (e.g., single streamtube model, multiple streamtube model, double multiple streamtube model) and vortex models (e.g., fixed-wake vortex model [25] and free-wake vortex model [17]).In streamtube models, the time-averaged forces on the blades are equated to the change in the streamwise momentum of the wind.Although, in principle, with vortex models, one can obtain the flow field downstream of a VAWT, the main focus of both models has been on the prediction of forces on the blades and overall rotor performance.
Simulation of VAWTs by solving the Navier-Stokes equations with numerical techniques was pioneered by Rajagopalan and Fanucci [26], who introduced the actuator swept-surface method in a two-dimensional space as a cylindrical porous shell on which the time-averaged blade forces continuously act on the fluid.This method was extended to a three-dimensional actuator swept-surface model for a curved-bladed VAWT in Rajagopalan et al. [27].In both studies, the flow was considered laminar.Fortunato et al. [28] applied the same turbine modeling and solved the 2D compressible laminar Euler equation with a finite difference method.Shen et al. [29] used a 2D Navier-Stokes solver, using the k − ω shear stress transport turbulence model, to obtain the flow field past a VAWT.They modeled the turbine blades as two rotating actuator surfaces, on which the body forces are distributed.In the above works, mainly power coefficients and forces acting on the blades are quantitatively compared with experimental data, while the wake flow field is only demonstrated qualitatively.Moreover, considering the highly turbulent wake of a VAWT, the application of a numerical turbulence-resolving technique can potentially improve the quality of the wake prediction.
The present work aims at simulating the turbulent wake of a VAWT in a large eddy simulation (LES) framework by solving the full filtered Navier-Stokes equations in three dimensions, using an actuator swept-surface and an actuator line method to model the turbine effects on the flow and two turbulence modeling techniques to parameterize the subgrid-scale (SGS) stresses.It is worth mentioning that it is the first time that LES has been applied and validated for the simulation of VAWT wakes by using either the actuator swept-surface model (ASSM) or the actuator line model (ALM) techniques.The LES framework and model formulations are presented in Section 2. Next, we simulate the wake behind a straight-bladed VAWT for one special case and compare the results with corresponding experimental data.The numerical and experimental setups are described in Section 3, and the results are presented in Section 4. A summary of the study is given in Section 5.

Large Eddy Simulation Framework
LES solves the filtered incompressible Navier-Stokes equations, which can be written in rotational form as: where the tilde represents a three-dimensional spatial filtering operation at scale ∆, u i is the filtered velocity in the i−direction (with i = 1, 2, 3 corresponding to streamwise (x), spanwise (y) and vertical (z) directions, respectively), p * = p ρ + 1 2 u i u i is the modified kinematic pressure where p is the filtered pressure, τ ij = u i u j − u i u j is the kinematic SGS stress, ν is the fluid kinematic viscosity, f i is a body force (per unit volume) used to model the effects of the turbine on the flow, F p is an imposed pressure gradient and ρ is the constant fluid density (following the Boussinesq approximation).For pointwise conservation of kinetic energy, the rotational form of the convective term is used (using the incompressible flow identity) [30][31][32].In the next subsections, we discuss the parameterization of SGS stresses (τ ij ) and the turbine-induced forces (f i ).

Turbulence Closure Models
Two closure models for parameterizing the SGS stresses are considered: the traditional Smagorinsky model and the modulated gradient model.

The Smagorinsky Model
The traditional Smagorinsky model [33] for the deviatoric part of the SGS stress is: where is the resolved strain-rate tensor, whose magnitude is | S|; ∆ is the filter width and C S is the Smagorinsky coefficient.As common practice in the application of the Smagorinsky model in the LES of wall-bounded flows, we use the ad hoc wall damping function proposed by Mason and Thomson [34] and commonly used for the calculation of the Smagorinsky coefficient near walls.Further details regarding the application of the Smagorinsky model in the present code (including the near-wall treatment) can be found, for instance, in Porté-Agel et al. [35] and Wan et al. [36].

The Modulated Gradient Model
The modulated gradient model (MGM), which was proposed by Lu and Porté-Agel [37], is an alternative to standard eddy-viscosity closure models.In the MGM, the SGS stress tensor is modeled as: If ∆ x , ∆ y and ∆ z are taken to be filter sizes in the x-, yand z-directions, then G ij is defined as: To evaluate the SGS kinetic energy, k sgs , the local equilibrium hypothesis is used, which assumes there is a local balance between SGS energy production, P , and dissipation, ε.These two terms are calculated as: and: where the coefficient, C ε , is assumed to have the theoretical value of 1.0 [38,39].Replacing τ ij with the expression of Equation ( 4) in Equation ( 6) and equating Equations ( 6) and (7), one can solve for k sgs : Replacing the above expression for k sgs in Equation (7) gives the following equation for the SGS dissipation, ε: To ensure numerical stability and that energy does not transfer locally from unresolved to resolved scales, we use the following clipping procedure: With Equation (10), the set of equations for the MGM model is closed, and τ ij is fully expressed in terms of resolved variables.The model has been successfully validated in the simulation of a neutral atmospheric boundary layer [37,40,41] and a boundary layer over a 2D block using the immersed boundary method [42].

Parameterization of VAWT-Induced Forces
In this study, two methods are used to model the turbine-induced forces: (a) the actuator swept-surface model (ASSM); and (b) the actuator line model (ALM).Both models are based on blade-element theory, and therefore, they do not require resolving the boundary-layer flow around the turbine blade surfaces, which greatly reduces the computational cost requirements.In the blade-element theory, the turbine blades are divided into N blade elements that are assumed to behave aerodynamically as two-dimensional airfoils.The aerodynamic forces are determined using the lift and drag characteristics of the airfoil type, as well as the local flow conditions.
Figure 1 shows a horizontal cross-section of a blade in Darrieus motion in an arbitrary angular position, θ.In this figure, Ω b is the rotor angular velocity, R is the radius of the rotor, u ∞ is the unperturbed wind velocity magnitude and V local is the local wind velocity vector at the position of the blade.Denoting the radial and tangential components of the local wind velocity in the inertial frame of reference as V local,n and V local,s , respectively, the local relative wind velocity vector with respect to the rotating blade can be written as , where e n and e s are unit vectors in the radial and tangential directions, respectively.As is common for VAWTs, the airfoil is considered to be symmetrical, and the chord-line is assumed to be normal to the rotor radius; thus the angle of attack can be written as: The lift and drag forces (see Figure 2) are calculated as: where  In the ASSM approach, the actuator swept surface is defined as the surface swept by the blades.For example, in the case of a straight-bladed VAWT, this actuator swept surface is an actuator cylinder.We distribute the time-averaged turbine-induced forces on this surface and integrate them over time and space.Figure 3 shows a horizontal cross-section of the actuator swept surface in the computational mesh, where the black points are grid points and the dashed lines show the control volumes surrounding them.

∆
For illustration reasons and for the sake of clarity, the size of the grid cells with respect to the diameter of the rotor is shown as being much bigger than in the simulations.The turbine forces are distributed among grid points, whose corresponding control volumes intersect the actuator swept surface.For instance, arc ab with an angle of ∆θ is located inside the control volume of grid point q in Figure 3.We take the midpoint, h, as a representative point of this arc and calculate V rel , C L and C D on this point.
The time fraction that each blade spends in arc ab is (∆θ/Ω b ) × (Ω b /2π), and considering the number of blades, B, the time-averaging weight for the calculated force on this point is B∆θ/2π.Finally, the term, f i /ρ, in the Navier-Stokes equations [Equation ( 2)] for a given grid point, q, is calculated as: where i and j are unit vectors in the streamwise and spanwise directions, respectively; and ∆x, ∆y and ∆z are grid sizes in the x-, yand z-directions, respectively.At each time-step, this calculation is done for all grid points whose control volumes intersect the actuator swept-surface, and the obtained forces are applied for each of them.It is worth mentioning that if one applies the ASSM approach to HAWTs, it would simply reduce to the conventional actuator disk model.

Actuator Line Model
In the ALM, the turbine forces are distributed vertically along lines that represent the blades of the wind turbine.In this model, which was introduced by Shen and Sørensen [43] and has already been implemented for horizontal-axis wind turbines [43][44][45][46][47], there is no need to time-average the forces, because the blades are tracked in each time-step, and the forces are applied only at the physical position of the blades.
First, each grid point whose corresponding control volume encompasses the blade line is identified, and then, the airfoil data and subsequent loading are determined by computing local angles of attack from the movement of the blades and the local flow field.For instance, in Figure 4, it can be seen that the blade is located at point h, which is within the control volume of grid point q; hence, theoretically, the term, f i /ρ, in the Navier-Stokes equations [Equation ( 2)] should only be assigned to this grid point with the following formula: By explicitly representing the individual blade motions, the ALM enables us to capture in more detail the complex and unsteady vortical structure of the near-wake flow.Consequently, it is likely to improve the near-wake predictions with respect to the ASSM.The application of the ALM seems particularly crucial for VAWTs, since, in contrast to HAWTs, a blade of a VAWT passes through the wakes generated by other blades, as well as its own wake; this fact makes the ALM approach attractive, because it allows us to study the blade-wake interaction phenomenon more realistically.

Numerical and Experimental Setups
In this section, the details of the numerical simulations, as well as the experimental setup, with which our results are compared, are presented.We have chosen the experimental study by Brochier et al. [18] to compare our simulation results with, because it contains data on mean velocity and turbulence intensity at several downstream positions in the wake of a model VAWT, even though it was conducted under uniform flow conditions.In fact, in previous VAWT wind tunnel experiments, the turbine has never been placed in a turbulent boundary layer flow, which is the regular flow regime, where full-scale turbines are located.It should be mentioned that more systematic and rigorous experimental studies on VAWT wakes, especially in a turbulent boundary layer flow, are still needed to deepen our understanding of the flow structure and behavior.

Experimental Setup
In 1986, Brochier et al. [18] conducted a water channel experiment with a straight-bladed Darrieus wind turbine under controlled-rotation conditions in the Institute de Méchanique et Statistique de la Turbulence (IMST) water channel.Velocity measurements were made using a laser-Doppler velocimeter (LDV), and quantitative profiles of mean velocity and turbulence intensity were provided.
The water channel was vertical, and the driving force of the flow was gravity.The length of the channel in the streamwise direction (vertical direction in the experiment) was 1.5 m, and the cross-section of the channel (normal to the streamwise direction) was a square of a side length of 20 cm (see Figures 5  and 6).The walls of the channel were transparent (made of plexiglass) to allow for LDV measurements and flow visualization.The model rotor was made up of two straight blades with a chord length of 2 cm.The blades are composed of a NACA 0018 airfoil.The length of the blades was 20 cm, and the diameter of the rotor was 12 cm.The center column of the turbine had a diameter of 1 cm.The measurements were obtained for two tip speed ratios of 2.14 and 3.85.The flow speed in the water channel was 15 cm/s, and the center of the turbine was located 13 chord lengths downstream of the nozzle plane through which the inflow enters the test section.

Numerical Setup
The LES code is a modified version of the code described by Albertson and Parlange [32], Porté-Agel et al. [35] and Porté-Agel et al. [47].A 3D structured mesh is employed, which has N x and N y nodes in the xand y-directions, respectively, and N z nodes in the z-direction.The mesh is staggered in the z-direction; this means the layers in which the vertical component of velocity (w) is stored are located halfway between the layers in which all the other main flow variables (u, v, p) are stored.The first w−nodes are located on the z = 0 plane, while the first uvp−nodes are located on the z = ∆z 2 plane.To compute the spatial derivatives, a Fourier-based pseudospectral scheme is used in the horizontal directions, and a second-order finite difference method is used in the vertical direction.The governing equations for conservation of momentum are integrated in time with the second-order Adams-Bashforth scheme.
The boundary conditions in the horizontal directions are mathematically periodic.For the bottom and top boundary conditions, the instantaneous surface shear stress is calculated using the Monin-Obukhov similarity theory [48] as a function of the local horizontal velocities at the nearest (to the wall) vertical grid points (z = ∆z 2 for the bottom wall and z = L z − ∆z 2 for the top wall) (see, for instance, [49][50][51]).In the simulations for this case, the same dimensions as in the experiment are used in the yand z-directions, while in the streamwise direction, the domain is taken longer than the channel length, because of the presence of a buffer zone (which will be described later).Thus, the domain dimensions are L x = 200 cm, L y = 20 cm (not including the width of the buffer zone for side-walls) and L z = 20 cm in the streamwise, spanwise and vertical directions, respectively.Regarding the computational mesh, the number of grid points in each of the three directions is N x = 240, N y = 30 (not including the 4 auxiliary grid points within the buffer zone for side-walls), and N z = 80.
Because of the periodic boundary conditions in the xand y-directions, it is implicitly assumed that there are similar turbines upstream, downstream and on both sides of the VAWT under consideration.To eliminate the effects of the upstream turbine wake on the flow and to ensure that we have a uniform inflow (as in the experiment), a buffer zone upstream of the turbine is employed to adjust the flow to an undisturbed uniform inflow condition.The use of this technique, i.e., using an inflow boundary condition in a direction in which the flow variables are discretized using Fourier series, has been shown to be successful in the works of Tseng et al. [52], Wan and Porté-Agel [53] and Wu and Porté-Agel [54].
In the y-direction, to account for the effect of the side walls, an immersed boundary technique is used.We have set the wall-normal velocity component (v) to zero at the location of the walls and have also modified the τ xy component of the stress tensor for the grid points, which are located on the side walls, to a value obtained from assuming a logarithmic velocity profile near the side walls.The computational domain is shown in Figure 7.In this figure, the location of Side-wall 2 can be assumed to be at the other end of the y-direction axis (because of the periodic boundary condition in y-direction), so that Side-wall 1 and Side-wall 2 represent the two side walls of the experiment, which are 20 cm apart.
As discussed in Subsection 2.2, for both VAWT models (ASSM and ALM), we need tabulated data for lift and drag coefficients.The conventional source for aerodynamic characteristics of symmetrical airfoils, used in VAWTs, is the tabulated data of Sheldahl and Klimas [55]; however, for low Reynolds numbers (which is the case for this experiment), these data show some anomalies for low angles of attack (i.e., the lift coefficient tends to be negative).This data anomaly has been reported by Shen et al. [29], Kumar et al. [56] and Hara et al. [57].As a solution, for low angles of attack, the data of Kumar et al. [56] is used, and for high angles of attack, the data of Sheldahl and Klimas [55] is used (for NACA 0018 airfoil).This approach to obtain the lift and drag coefficients has been employed by Hara et al. [57], as well.To account for the dynamic stall phenomenon and to consider its effect on C L and C D curves, the modified MIT dynamic stall model is employed [58].It is worth mentioning that the the dynamic stall phenomenon has a significant effect on the values of the airfoil lift and drag coefficients and, hence, on the turbine performance [1,29,59].The turbine center mast is also modeled as a cylindrical body, which applies a drag force on the fluid.The drag coefficient for the mast is considered to be equal to 1.0.In our simulations, the tip speed ratio (defined as the ratio between the velocity of the blades and the incoming flow velocity) of the turbine is taken to be 3.85 to match that of the experiment.The inflow of the simulation is a uniform flow with a velocity of 15 cm/s and a turbulent intensity of 2.5%, as in the experiment.

Results and Discussion
In this section, the results of our simulations for the above-described case are presented for five different combinations of the turbine and SGS models.Three simulations are performed for the ALM (with the MGM and the traditional Smagorinsky model with C S = 0.1 and 0.2), and two are done for the ASSM (with the MGM and the traditional Smagorinsky model with C S = 0.1).
Figures 8 and 9 show time-averaged streamwise velocity contours on a horizontal plane (top view) at the equator height of the turbine for the actuator swept-surface and actuator line models, respectively.The black circle in each figure represents the flight path of the turbine blades, which rotate in the counter clockwise direction.It is obvious from the figures that the turbine extracts momentum from the flow, resulting in a low velocity region in its wake.As can be seen in the contours, when we go further downstream of the turbine, the wake recovers by entraining the neighboring high velocity fluid.In these figures, the blockage effect is evident near the channel side-walls (two ends of the y-axis), where the streamwise velocity magnitude is substantially larger than the incoming flow velocity.This is due to the narrow channel in which the experiment was performed.It should be noted that the presence of the blockage effect affects the recovery rate of the wake, due to the higher velocity gradients at the edge of the wake with respect to unconfined VAWTs.In these velocity maps, the effect of the turbine wake can be seen at downstream distances as large as 20 rotor radii.In all the following figures the streamwise (X) and spanwise (Y ) distances are non-dimensionalized by the rotor radius, R, and the velocities are normalized by the incoming uniform flow velocity, U 0 .As can be observed in the figures, the wake of a VAWT is not symmetric in the spanwise direction (y-direction).This lack of symmetry, which is a typical behavior of VAWT wakes, was also reported by Battisti et al. [20], Simao Ferreira et al. [60] and Brochier et al. [18].We can see that in the region where the streamwise component of the blade velocity is in the opposite direction of the inflow velocity, i.e., upcoming blade region, (positive y in the figure), the wake is stronger than in the region where these two velocities are in the same direction, i.e., retreating blade region, (negative y in the figure).The primary reason behind this observation can be found by inspecting the velocity triangles at two typical blade positions in these two regions, as shown in Figure 10.In Figure 10a, the black circle is the flight path of the blades, on which four typical blade positions are highlighted.From the point of view of an observer who sits on the blades, the blade velocity is a constant vector (with a magnitude of RΩ b ) and the induced flow velocity at the point of the blade is a vector that rotates with the same average frequency as Ω b .The dotted circle in Figure 10b is the locus of the endpoint of the induced velocity vector, and the four points that correspond to the highlighted blade positions in Figure 10a are illustrated on this circle, as well.It is clear from the velocity triangles of points a and b that the relative velocity at point a is larger than the one at point b.As the blade forces are proportional to the square of the relative velocity at the blades, this explains why the velocity deficit is stronger at points a and d as compared to points b and c.It is noteworthy that in this analysis, for the sake of simplicity, the induced velocities at all blade positions are considered to be the same; however, it is evident that the induced velocities at points d and c are smaller than the ones at points a and b.This is because the downwind half of the rotor is located in the wake of the upwind half.To obtain a more quantitative understanding of the turbine wake, mean velocity profiles in the spanwise direction and at the equator height of the rotor are plotted in Figure 11 for five downstream locations: x/R = 2.50, 3.33, 5.0, 8.33, 11.67.In these plots, the five combinations of SGS and VAWT models are compared together and also with experimental data.As can be seen in the figure, although both turbine models are able to predict the evolution of the wake in a fairly acceptable manner, the agreement between the ALM results and the experimental data is slightly superior to that of the ASSM results.One important feature that can be observed in this figure (and also in Figures 8 and 9) is the different ability of each model in predicting the wake recovery rate.Clearly, the wake recovers more slowly in the simulations run with the ASSM than in the ones with the ALM.This relatively slow wake recovery is also observed for the Smagorinsky model with C S = 0.2.Since the wake recovery rate is closely correlated to the turbulence intensity level in the flow (e.g., Wu and Porté-Agel [61] for horizontal-axis wind turbines), it would be insightful to look at how the VAWT generates turbulence in the flow.Figures 12 and 13 show contours of the mean streamwise turbulence intensity (TI) in a horizontal plane at the equator height of the rotor for the ASSM and the ALM, respectively.As can be understood readily from these contours, the difference between the ASSM and the ALM is even more pronounced considering the distribution of turbulence intensity.In these figures, one should pay attention to the trends in turbulence intensity distribution in the spanwise and streamwise directions.In the spanwise direction, all of the models qualitatively predict two regions of relatively high turbulence intensity in the vicinity of the wake edges.More details on the spanwise variation of TI will be discussed later.While in the spanwise direction, the same trends exist, in the streamwise direction, the differences between different models are more obvious.As can be seen in the contours, for the ASSM model, we have very low levels of TI in the turbine region and also in the near wake; then, TI gradually increases, until a point where it reaches a maximum (x/R around 10), and then, it starts to decrease.However, for the ALM, we observe a very high TI in the turbine region and in the near wake, and it decreases as we go downstream.Of these two trends, the one obtained by the ALM is qualitatively in better agreement with the experiments of Brochier et al. [18].Looking at instantaneous velocity fields (Figure 14) can shed more light on the dissimilar effects of the two turbine models on the near wake flow.In Figure 14b (the ALM simulation), the wavy and wiggling edges of the wake near the turbine are clearly distinguishable from the rather straight and smooth wake edges in the same region in the ASSM simulation (Figure 14a).One can attribute this difference to the fact that the ALM can better capture the vorticity shed from the blades in the near wake by resolving the motion of individual blades.Another interesting observation that can be made in Figure 13 is the existence of four locations on the flight path of the blades around which the TI is maximum.This is likely due to the fact that, in those locations, the magnitude of the forces induced by the flow on the blades is locally maximum.In Figure 15, where the spanwise variation of TI is shown for five downstream locations (x/R = 2.50, 3.33, 5.0, 8.33, 11.67), the above-discussed streamwise trends can be seen more clearly for the ASSM (blue lines) and the ALM (black lines).Also seen in this figure is the fact that in the spanwise direction the TI profile has two maxima, which are located on the edges of the wake.This is due to the unsteady shed vorticity from the blades and is reported in the works of Fujisawa and Shibuya [21], Battisti et al. [20] and Brochier et al. [18].Moreover, we see different behavior of the flow turbulence in the spanwise direction for different turbine models.The ASSM overestimates the difference of TI between the retreating-blade half of the domain (negative y) and the upcoming-blade half of the domain (positive y); in contrast, the TI profiles obtained from the ALM show a better agreement with the measurements.This can be attributed in part to the inability of the ASSM to capture the unsteady-periodic nature of the flow, especially in the near wake region.This is, in contrast to the ALM, that, by resolving the unsteadiness of the blade forces, which improves the prediction of the turbulent behavior of the flow, particularly in the retreating-blade half of the domain.
At this point, we are able to draw conclusions regarding the links between the mean velocity field and the turbulence predictions of different models.First of all, as discussed earlier, the fact that the ASSM predicts lower levels of TI in the near wake leads to slower wake recovery.This can be seen clearly in Figure 11e, where the ASSM simulations (blue lines) overestimate the velocity deficit at x/R = 11.67 with respect to the ALM simulations (black lines) and also the measurements.Regarding the SGS models, we can see that, as expected, the high SGS dissipation rate forced by the Smagorinsky model with C S = 0.2 leads to lower resolved turbulence intensity predictions (Figure 15, green line) and, consequently, to a slower wake recovery (Figure 11e, green line).
By observing the presented results, one can conclude that the ALM, coupled with the MGM or the Smagorinsky model with C S = 0.1, gives the best agreement with the measurements.However, it should be noted that the MGM has a clear advantage over the Smagorinsky model in that it does not need tuning for its model coefficient [41].In contrast, even though the value of C S = 0.1 gives the best agreement of the Smagorinsky model results with experimental data for this specific problem, the optimum value could be different for other experimental or field conditions.Specifically, this advantage of the MGM is expected to be even more decisive in the case of VAWTs placed in turbulent boundary layer flows, for which the MGM has been shown to provide better predictions (without turbines) than the standard Smagorinsky model [37,40,41].Another advantage of the MGM over the Smagorinsky model (in fact, all Smagorinsky-based models) is the ability of the MGM to calculate the unresolved part of the TI in addition to its resolved part.This is due to the fact that in the MGM, the complete SGS stress tensor is modeled, which is in contrast to Smagorinsky-based models for which only the deviatoric part of the SGS stress tensor is parameterized.In our case, the unresolved part of the TI was about 5% of the total TI throughout the domain.To show the grid resolution sensitivity of the computational results, we performed simulations for the ALM-MGM and ASSM-MGM model combinations with three different grid resolutions (the resolution is varied in the horizontal directions, which is critical in capturing the blade motions).It was observed that, for the range of resolutions considered here, resolution had only a small effect on the results.Figure 16 shows spanwise plots of mean streamwise velocity and turbulence intensity in the middle downstream position (x/R = 5.0) for the three grid resolutions.

Summary
In this study, large eddy simulation (LES), coupled with a turbine model, is used to investigate the 3D unsteady turbulent wake flow behind a vertical-axis wind turbine (VAWT).To model the SGS momentum fluxes, two SGS models are used, namely, the traditional Smagorinsky model and the modulated gradient model (MGM) [37].To parameterize the effects of a VAWT on the flow, two VAWT models are developed: (a) the actuator swept-surface model (ASSM), in which the time-averaged turbine-induced forces are distributed on a surface swept by the turbine blades, i.e., the actuator swept-surface; and (b) the actuator line model (ALM), in which the instantaneous blade forces are only spatially distributed on lines representing the blades, i.e., the actuator lines.In both models, blade-element theory is used to calculate the lift and drag forces acting on the blades.This is the first time that LES has been applied and validated for the simulation of VAWT wakes by using either the ASSM or the ALM technique.The results are compared with flow field measurements in the wake of a model straight-bladed VAWT, carried out in the Institute de Méchanique et Statistique de la Turbulence (IMST) water channel.
Different combinations of VAWT and SGS models are tested and compared against each other, and a fairly good agreement with the experimental data is observed.It is found out that, in general, the results obtained with the ALM agree better with the measurements and also with the physics of the flow than the ones from the ASSM.The inherent advantage of the ALM in resolving the periodic location of the blade forces allows one to capture the unsteady-periodic nature of the wake and, thus, the blade-wake interaction, which is a characteristic of the aerodynamics of a VAWT rotor.The slight superiority of the ALM over the ASSM was shown quantitatively in the better prediction of the turbulence intensity and, consequently, of the wake recovery, which is very important in designing and optimizing potential VAWT arrays.However, it is worth bearing in mind that the ASSM could potentially be the preferable turbine model in simulations of VAWT wind farms, due to its lower grid resolution requirement.Regarding the SGS models, although both tested models showed reasonable turbulence modeling capabilities, it should be noted that the traditional Smagorinsky model needs its model coefficient to be optimized, depending on the local flow characteristics, while the MGM proved that it can, with its theoretically-determined model coefficient, capture the mean flow and turbulent characteristics of a VAWT wake.This advantage of the MGM is expected to be even more decisive in the case of VAWTs placed in turbulent boundary layer flows, where the MGM has previously been shown to have a superior prediction ability (without turbines) compared with the standard Smagorinsky model [37,40,41].
Potential sources of the observed discrepancies between simulation results and measurement data can be summarized as: (a) inherent limitations of the proposed numerical model; (b) inherent measurement errors in the experiment; and (c) inaccuracies in the tabulated airfoil data.It is worth noting that, in effect, the final results are considerably dependent on the lift and drag coefficient curves as a function of the angle of attack.
As mentioned by Paraschivoiu [1], the most important shortcoming that VAWTs are suffering from, at present, is that they have not benefited from research and development as much as HAWTs.Therefore, more rigorous and systematic studies on different technical aspects of VAWTs are still needed.For the aerodynamic aspect, high-resolution wind tunnel measurements of the flow field in the wake of VAWTs placed in a boundary layer flow would be highly valuable.Furthermore, for the numerical simulation of the wake, the authors believe that, if one lesson has to be learned from the development of HAWTs, it is that the most feasible and practical way of taking into account the effects of the turbine on the flow is through using immersed-body-force approaches, such as an actuator swept surface, actuator line or rotating actuator surface.With these techniques, the simulation, optimization and design of VAWT wind farms will be an attainable and realizable objective in the near future.
Re c ) and C D = C D (α, Re c ) are the lift and drag coefficients obtained from tabulated airfoil data, respectively, Re c is the Reynolds number based on relative velocity and chord length, e L and e D are unit vectors in the lift and drag directions, ∆l is the vertical length of the two-dimensional airfoil element and c is the chord length.

Figure 3 .
Figure 3. Schematic of a cross-sectional view of the actuator swept-surface in the computational mesh.

Figure 4 .
Figure 4. Schematic of a cross-sectional view of the actuator line in the computational mesh.

Figure 5 .
Figure 5. Schematic of y-z view (normal to streamwise direction) of the experimental setup.

Figure 6 .
Figure 6.Schematic of x-z view of the experimental setup.

Figure 7 .
Figure 7.  x-y view (normal to the blades) of the computational domain.

Figure 8 .Figure 9 .
Figure 8. Contours of the normalized mean streamwise velocity, u/U 0 , in a horizontal plane at the equator height of the rotor for the actuator swept-surface model (ASSM): (a) Smagorinsky model, Cs = 0.1; (b) modulated gradient model (MGM).

Figure 10 .
Figure 10.(a) A cross-sectional view of the flight path of the blades (black circle) and four azimuthal blade positions (colored circles); (b) the locus of the endpoint of relative velocity vector (dotted circle) and the corresponding velocity triangles.

Figure 12 .Figure 13 .Figure 14 .
Figure 12.Contours of the normalized mean streamwise turbulence intensity (resolved part), σ u /U 0 , in a horizontal plane at the equator height of the rotor for the ASSM: (a) Smagorinsky model, Cs = 0.1; (b) MGM.