Impact of Subgrid-Scale Modeling in Actuator-Line Based Large-Eddy Simulation of Vertical-Axis Wind Turbine Wakes

: A large-eddy simulation (LES) study of vertical-axis wind turbine wakes under uniform inﬂow conditions is performed. Emphasis is placed on exploring the effects of subgrid-scale (SGS) modeling on turbine loading as well as on the formation and development of the wind turbine wake. In this regard, the validated LES framework coupled with an actuator-line parametrization is employed. Three different SGS models are considered: the standard Smagorinsky model, the Lagrangian scale-dependent dynamic (LSDD) model, and the anisotropic minimum dissipation (AMD) model. The results show that the SGS model has a negligible effect on the mean aerodynamic loads acting on the blades. However, the structure of the wake, including the mean velocity and turbulence statistics, is signiﬁcantly affected by the SGS closure. In particular, the standard Smagorisnky model with its theoretical model coefﬁcient (i.e., C S ∼ 0.16) postpones the transition of the wake to turbulence and yields a higher velocity variance in the turbulent region compared to the LSDD and AMD models. This observation is elaborated in more detail by analyzing the resolved-scale turbulent kinetic energy budget inside the wake. It is also shown that, unlike the standard Smagorinsky model, which requires detailed calibrations of the model coefﬁcient, the AMD can yield predictions similar to the LSDD model for the mean and turbulence characteristics of the wake without any tuning.


Introduction
Power harvesting from atmospheric turbulent flows can be achieved through both horizontal-axis and vertical-axis wind turbines. Horizontal-axis wind turbines (HAWTs) are dominating the market, but vertical-axis wind turbines (VAWTs) have attracted renewed consideration in recent years due to their potential to provide a higher power output per unit land area [1]. This is, in part, because unlike for HAWTs, it is possible to increase the swept area of VAWTs independently of their footprints [2].
In addition to laboratory experiments [3][4][5][6][7][8][9][10][11] and field-scale measurements [2,12,13], numerical studies have also been performed in recent years to investigate the wake aerodynamics behind VAWTs [14][15][16][17][18][19]. Numerical simulations are a cost-effective approach that can provide detailed information about the fluid mechanics of VAWTs over a wide range of temporal and spatial scales. Specifically, large-eddy simulations (LES) coupled with a turbine model, such as the actuator-line technique, has become the leading numerical technique to simulate velocity and turbulence fields through wind turbines (see reviews by [20][21][22]). In the LES approach, turbulent structures larger than a certain filter width are resolved and the influence of the unresolved scales on the resolved motions are represented by a subgrid-scale (SGS) model. Hence, the accuracy of the LES can be affected by the performance of the SGS parametrization, especially in the anisotropic turbulence. It is of particular interest for wind energy applications to better understand and quantify how the SGS models affect the aerodynamic loads acting on the blades and also the formation and recovery of turbine wakes. This effect has been well studied and documented in the LES of HAWT wakes [23][24][25][26][27]. However, in the case of VAWTs, the role of SGS modeling on both the turbine loading and wake-flow structures has been studied less in prior investigations.
The present work aims to explore the effect of SGS modeling in LES on the turbine loading and flow characteristics of the wake downwind of a full-scale VAWT under uniform inflow conditions. The main reason for considering a uniform inflow condition in this study is that, in the absence of mean wind shear and turbulence in the incoming flow, the SGS model mainly controls the transition of the wake to turbulence [24,27]. Hence, this case will better assess the capability of different SGS models in capturing the point of transition and breakdown location of the wake. In this regard, the validated LES framework coupled with an actuator-line parametrization for the wind turbine is used [18]. A detailed description of the LES framework and the numerical setup are given in Section 2. The effects of SGS models on the aerodynamic loads as well as the spatial distribution of the mean and turbulence fields inside the wake are provided and discussed in Section 3. A detailed analysis of the resolved-scale turbulent kinetic energy (TKE) budget in the wake is presented in Section 4. Finally, concluding remarks are included in Section 5.

LES Governing Equations
The LES framework utilized in this study is based on the filtered incompressible Navier-Stokes equations and the filtered continuity equation, where the tilde indicates spatial filtering at the ∆ scale. ρ is the fluid density.ũ i is instantaneous filtered velocity where i = 1, 2, 3 corresponds to the streamwise (x), spanwise (y) and wall-normal (z) directions, respectively. t is time,p is the filtered kinematic pressure, τ ij = u i u j −ũ iũj is the SGS stress tensor, and f i is an external body force corresponding to the effect of wind turbine on the flow. The viscous stress tensor based on filtered variables is defined as σ ij = 2νS ij where ν is the fluid molecular viscosity, andS ij = (∂ iũj + ∂ jũi )/2 is the resolved rate-of-strain tensor.

Subgrid-Scale Parametrization
In LES, the dynamics of the large-scale motions are computed explicitly, and the influence of the small-scale (or unresolved) structures on the large scales is parametrized using an SGS model [28]. An important category of SGS closures consists of eddy viscosity models which take into account the effect of small-scale eddies by locally increasing the molecular viscosity by means of an eddy viscosity. In this approach, the SGS stress tensor is represented as where ν e is the eddy viscosity. In this study, different SGS models were considered, all of them within the class of eddy-viscosity closures: the standard Smagorinsky model, the Lagrangian scale-dependent dynamic model, and the anisotropic minimum dissipation model. A short explanation of these models is provided in the following section.

Standard Smagorinsky Model
Using the mixing length approximation is one possible way to determine the eddy viscosity. This approximation yields the well-known Smagorinsky model [29]. In this approach, the eddy viscosity is formulated as where C S is the Smagorinsky coefficient and |S| = 2S ijSij . In the standard Smagorinsky model (denoted SMAG below), the model coefficient is assumed to be constant, independent of space and time. For the case of homogeneous isotropic turbulence with ∆ within the inertial subrange, an analysis performed by Lilly [30] using the spectral cut-off filter yielded C S ∼ 0.16 . However, assuming a constant value for the model coefficient raises important concerns in the simulation of complex flows. In particular, for cases involving inhomogeneous and transitional flows, the general conclusion from a posteriori testing is that the standard Smagorinsky model is too dissipative; thus, the model coefficient should generally be decreased [28,[31][32][33][34]. It follows that, in practice, and in order to yield acceptable results, the standard nondynamic approach requires detailed calibrations of the model coefficient [35,36]. Here, to assess the effect of Smagorinsky coefficient on wake flow structures, two values for the model coefficient were considered. The first one was the theoretical value proposed by Lilly [30] for isotropic turbulence (i.e., C S ∼ 0.16) and the second one was the value used by Moin and Kim [32] for inhomogeneous channel flows (i.e., C S ∼ 0.065).

Lagrangian Scale-Dependent Dynamic Smagorinsky Model
The dynamic procedure introduced by Germano et al. [37] provides a methodology to determine the appropriate local value of the Smagorinsky coefficient without any ad hoc tuning. In this approach, the Smagorinsky coefficient, which is space and time dependent, is computed by comparing the eddy dissipation at different scales. The dynamic model has been shown to improve SGS dissipation characteristics compared to the SMAG model [38][39][40]. In addition, the model appropriately switches off and does not yields extra eddy dissipation in laminar and transitional flow regions. However, an important assumption in the standard dynamic model is that the model coefficient is scale invariant. This assumption was later relaxed through the development of a scale-dependent dynamic model [35]. In this study, the Lagrangian scale-dependent dynamic (LSDD) model is employed where the model coefficient is computed dynamically and averaged along the fluid particle path [41]. Note that the computational complexity of the LSDD is much higher than that of the SMAG model. This results from performing the Lagrangian averaging procedure as well as the application of two test-filter operators. Specifically, as documented in [42], the application of the LSDD model increases the computational cost by about 35-40% over the SMAG model. More details about the LSDD closure can be found in [42].

Anisotropic Minimum Dissipation Model
Minimum dissipation models provide an alternative approach to specify the SGS stress tensor. They give the minimum amount of eddy dissipation needed to dissipate the energy of the subgrid scales. Verstappen [43] developed the first minimum dissipation model (called the QR model where Q and R represent the invariants of the resolved strain-rate tensor) for the eddy viscosity. The QR model yielded accurate results in simulations of turbulent flows with isotropic grids. Later, Rozema et al. [44] introduced the anisotropic minimum dissipation (AMD) model in which the properties of the QR model are generalized to anisotropic grids. The extension of the AMD closure to modeling the buoyancy-driven turbulent flows was recently done by Abkar et al. [45,46]. Among many desirable theoretical and practical properties of the AMD model, it appropriately switches off in laminar or transitional flow regions, and it is more cost-effective than the dynamic Smagorinsky approach. The AMD eddy-viscosity model is given by is the scaled gradient operator, in which C µ is a modified Poincaré constant, and ∆ µ is the grid size (Greek subscripts are excluded from the summation convention). Note that the modified Poincaré constant only depends on the order of accuracy of the discrete derivative operator in each direction, and its magnitude is independent of the size of the grid.
In particular, the Poincaré constant is 1/ √ 3 for second-order accurate schemes, and is equal to 1/ √ 12 for spectral methods [44,47]. The AMD model showed robust performances in simulations of temporal mixing layers, decaying grid turbulence, and turbulent channel flow [44], in simulations of thermally stratified atmospheric boundary layers [45,46], and in simulations of VAWT wakes [18]. More details on the derivation of the AMD model for the SGS stress tensor can be found in [44,46].

Wind Turbine Parametrization
In this study, the actuator-line model was used to parametrize the turbine-induced forces. Through this approach, the drag and lift forces acting on the blades were computed using the blade element theory and applied along the lines which represent the blades of the turbine. In Figure 1, a horizontal cross-section of the blades is shown. Different velocities, angles, and forces are shown in this figure.
Ω is the turbine angular velocity, and R is the radius of the rotor. V local is the local velocity of the incident flow at the blade. α is the angle of attack, and V rel is the local velocity relative to the rotating blade. The turbine-induced force per unit vertical length of the blade is provided by where C D = C D (α, Re c ) and C L = C L (α, Re c ) are the drag and lift coefficients, respectively; Re c is the Reynolds number based on the chord length, c, and relative velocity, V rel ; e D and e L denote the unit vectors in the directions of the drag and lift forces, respectively. The effect of dynamic stall on the drag and lift coefficients is taken into account using the modified MIT dynamic stall model [48]. This method has been successfully tested in recent studies focusing on modeling of VAWT wakes [10,14,15,18]. A tip-loss correction similar to the one proposed by Glauert [49] was also used in the blade element analysis. During the simulation, the local velocity V local at the blades was known. As a result, the local velocity relative to the rotating blades, V rel , and the angle of attack, α, could be calculated. The resulting forces were obtained using Equation (5). Note that, for numerical stability reasons, the applied blade forces were distributed smoothly over the actuator lines through the convolution of the computed force, f, and a regularization kernel, η [50], as where d is the distance between the actuator-line and computational grid points, and is a parameter that controls how the forces are distributed onto the grids. The value of is typically in the range of 1-3 grid sizes [51]. In this study, it was set to 1.5 times the grid size. The effect of the value on the simulation results was further assessed and is presented in Appendix A. Note that the present LES framework using the actuator-line technique has been validated against the laboratory data in the wake of a straight-bladed VAWT placed in the water channel [18].

Numerical Setup
The research code presented here has been developed and used in a series of wind energy studies (e.g., [52][53][54][55][56][57][58]). The numerical discretization follows the approach used by Moeng [59] and Albertson and Parlange [60] in which the computational domain is divided uniformly into N x , N y , and N z grid points in the x, y, and z directions, respectively. The governing equations (i.e., Equation (1)) are solved numerically using a pseudo-spectral discretization in the horizontal directions and a second-order finite difference scheme in the wall-normal direction. Aliasing errors are corrected in the nonlinear term using the 3/2 rule [61]. For time advancement, a second-order accurate Adams-Bashforth method is employed. The boundary conditions are periodic in the horizontal directions. Zero stress and zero vertical velocity conditions are imposed in the top and bottom boundaries.
In the simulations, a 200 kW VAWT with a rotor diameter (D) of 26 m was immersed in the flow. This particular turbine, which was designed and manufactured in collaboration with Uppsala University [62], is a straight bladed, Darrieus-type turbine. The rotor includes three straight blades with a span (H) of 24 m. A standard NACA 0018 airfoil with a chord length of about 0.75 m was used for the blades [63]. Below the rated wind speed (∼12 m/s), the tip-speed ratio of the turbine was set to 3.8. The drag and lift coefficients of the blades were obtained from the tabulated airfoil data provided in [64]. The averaged Reynolds number based on the chord length and relative velocity was about 1.52 × 10 6 for a tip-speed ratio of 3.8.
The domain sizes in the x, y, and z directions were set to L x = 30D, L y = 6D, and L z = 6D, respectively, and the simulations were performed with a grid resolution of 480 × 96 × 192. Note that the grid resolution sensitivity of the numerical results was tested and is provided in the Appendix A. The turbine was placed at three rotor diameters downstream from the inflow in the center of the y − z plane. The incoming wind to the turbine was uniform with a velocity of U 0 = 8 m/s. To adjust the flow from the downstream wake state to an undisturbed inflow condition, a fringe zone upstream of the turbine was used. In this method, the downstream wake flow does not affect the inflow condition due to the imposed periodic boundary conditions [65][66][67][68]. The simulations were initially run for two flow throughs to guarantee that quasi-steady conditions were achieved, and then the quantities of interest were averaged over three flow throughs. Figure 2 shows the time-averaged values of the angle of attack and the tangential force coefficient versus the azimuthal angle, respectively. The results obtained from the blade element momentum (BEM) theory are also shown for comparison. Note that the BEM theory equations were derived from momentum and blade element theories and are solved iteratively [69]. As can be seen in these figures, there was good agreement between the LES simulations and the BEM theory in the first half of the blade revolution. However, there was a discrepancy between the LES and the BEM theory in the second half of the blade revolution. This is mainly because in the second half of the revolution, the blades are in the wake generated by the other blades in the upstream. Since the reference velocity used in the BEM theory is independent of the blade locations, this approach is not able to correctly predict the aerodynamic loads in the second half of the blade revolution. These figures also reveal that the SGS model negligibly affects the mean aerodynamic loads acting on the blades and the predictions of all SGS models lay on top of each other.  Note that the wake generated behind the turbine gradually recovers with distance as faster moving fluid is entrained in the wake. The asymmetry of the wake in the horizontal plane can also be observed in this figure which is the inherent characteristic of VAWT wakes and is mainly related to the rotation of the turbine. Similar behavior has been reported in previous numerical and experimental studies of the wake of VAWTs [3,5,8,17,[70][71][72][73]. The vertical and lateral profiles of the mean streamwise velocity at the turbine equator are also plotted in Figure 4 for five different downstream locations at x/D = 1, 3, 7, 15 and 24. It can be visually acknowledged that the flow fields predicted by all the SGS models are similar in the region close to the turbine. The differences between the profiles are observed at locations where the turbulence is trigged and the wake starts breaking down. In particular, the SMAG model with the theoretical value for the Smagorisnky coefficient (C S ∼ 0.16) postpones the breakdown location of the wake. This is most likely because the model is too dissipative and damps the small-scale turbulence that causes the instability of the wake [26,27]. It was also found that the AMD model can provide comparable results to the LSDD model for predicting the wake velocity. Note that the computational complexity of AMD is much lower than LSDD, and it is comparable to the SMAG model. Additionally, it was realized that decreasing the model coefficient in the SMAG approach to C S ∼ 0.065 improves the mean velocity predictions in the vicinity of the breakdown location of the wake.   Figure 5 shows the contour plots of normalized streamwise velocity variance u u /U 2 0 in the central vertical and horizontal planes through the turbine. Here, the overbar denotes temporal averaging, and u =ũ −ū. The lateral and vertical profiles of this quantity through the center of the turbine are also plotted in Figure 6. It can be observed that in the near wake region in which the turbulence is not triggered yet, the velocity variance is almost zero. At some distance downwind from the turbine where the turbulence is triggered, an apparent enhancement of turbulent fluctuations is observed. Specifically, a higher turbulence intensity is observed at the top and bottom edges and close to the two side-tip positions of the turbine, which is mainly caused by the strong shear at these locations [15]. It can also be observed that, in the horizontal plane, the turbulence level is stronger in the upcoming blade region which is connected to stronger shear in that area compared to the retreating blade region. These figures also investigate the effect of the SGS model on the wake formation and breakdown. It can be understood from these figures that the location at which the transition of the wake to turbulence occurs is affected by the SGS closure. Consistent with the previous results for the mean velocity profiles, the SMAG model with C S ∼ 0.16 delays the point of transition which is due to the fact that the model is highly dissipative [27]. Note that unlike the mean velocity profiles, in which the predictions of different SGS models are similar in the far wake region, the discrepancies in the velocity variance profiles are larger and clearly visible at locations further downstream of the wake breakdown. It is also interesting to observe that despite the fact that the SMAG model with C S ∼ 0.16 is too dissipative, the turbulence level in the far wake is higher for that model than for the other cases which are less dissipative. This observation is investigated in more detail in the next section by performing the resolved-scale TKE budget analysis. Again, it can be noticed that the results obtained from the AMD model are comparable with the ones obtained from the LSDD model. Also, similar to the presented results for the mean velocity distribution, reducing the model coefficient in the SMAG model can improve the prediction of the turbulence properties in the far-wake region.

Resolved-Scale Turbulent Kinetic Energy Budget
As shown in Figure 5, the turbulence level in the far-wake is larger for the SMAG model with C S ∼ 0.16 than for the other cases. To better understand the role of the SGS closure on the energy exchanges between the resolved and the SGS parts, a detailed analysis of the resolved-scale TKE budget in the wake was carried out. The resolved-scale TKE budget equation can be expressed as [74] ∂ tk wherek denotes the mean resolved-scale TKE. T t , T p , T ν , and T sgs represent the transport of the resolved-scale TKE with velocity, pressure, viscous stress fluctuations, and SGS stress fluctuations, respectively. P s is the shear production describing the conversion of the mean kinetic energy to the resolved-scale TKE, P t is the rate of the work against the turbine-induced forces, ν is the viscous dissipation rate of the resolved-scale TKE, and sgs is the SGS dissipation rate and represents the energy transfer from the resolved scales to the subgrid scales. In this study, the focus was placed only on the last term on the right-hand side of the above-mentioned equation to investigate the contribution of the SGS closure on the energy exchange between the resolved and subgrid scales. Figure 7 displays the two-dimensional contours of the SGS dissipation rate obtained from the application of different SGS closures. It was observed that the SGS dissipation rate inside the wake is strongly affected by the SGS model, and in the particular cases considered here, it was lower for the SMAG model with C S ∼ 0.16 than for the other models. A possible explanation for this observation is that the energy transfer between the resolved scales is flawed in the SMAG model with C S ∼ 0.16, and that leads to the accumulation of kinetic energy in the resolved part and consequently, to a higher velocity variance in the turbulent region of the wake (shown in Figure 5). It is worth mentioning here that similar behavior was reported recently in the study by Martinez-Tossas et al. [27] in the simulation of HAWT wakes under uniform inflow. To justify this trend, they argued that C S ∼ 0.16 enhances the damping of the smallest resolved turbulent eddies which, in turn, disrupts the natural energy cascade and leads to the augmentation of TKE in the resolved part. The results of the TKE distribution in the wake are displayed in Figure 8.

Conclusions
The validated LES framework coupled with the actuator-line method was employed to investigate the impact of SGS modeling on the turbine loading as well as the formation and recovery of the wake behind a full-scale VAWT in uniform inflow. Three different SGS models were considered: the standard Smagorinsky (SMAG) model with two different values for the model coefficient, the Lagrangian scale-dependent dynamic (LSDD) model, and the anisotropic minimum dissipation (AMD) model.
The simulation results show that the influence of SGS models on the mean aerodynamic loads on the blades is negligible. However, the SGS model has a significant impact on the distribution of the wake velocity and turbulence statistics downwind of the turbine. In particular, the SMAG model with its theoretical value (i.e., C S ∼ 0.16) postpones the transition of the wake to turbulence to a distance further downstream and produces a higher turbulence intensity in that region compared to the LSDD and AMD models. A detailed analysis of the resolved-scale TKE budget in the wake region also revealed that the SMAG model with C S ∼ 0.16 overdissipates the energy of the smallest resolved scales. This, in turn, reduces the the natural energy transfer from the resolved to subgrid scales and consequently, leads to the accumulation of TKE in the resolved part [27]. It was also shown that calibrating the model coefficient in the SMAG model (i.e., C S ∼ 0.065 in this study) can improve the prediction of the turbulence statistics as well as the breakdown location of the wake. However, the AMD model can yield acceptable results similar to the LSDD model for the mean and turbulence characteristics in the wake without any tuning and with low computational complexity.
Funding: This research received no external funding.
Acknowledgments: This paper benefited from comments by Prof. Parviz Moin. The computational resources used in this study were provided by the Center for Turbulence Research (CTR) at Stanford University.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Effects of the Smearing Parameter and Grid Resolution
As mentioned in Section 2.3, for the sake of numerical stability, the applied blade forces were smoothly distributed through the actuator lines. In the application of the actuator-line model, a value in the order of 1 − 3 grid sizes is recommended for the smearing parameter ( ) shown in Equation (6) [51,75]. In order to study the sensitivity of the numerical results to the smearing parameter, numerical experiments were performed with two different values both within the recommended range: = 1.5∆x and = 2.5∆x. The lateral profiles of the normalized mean velocity and streamwise velocity variance obtained from the LSDD model with two different values are shown in Figure A1. As can be noticed from this figure, for the particular cases considered in this study, the profiles of the mean velocity and velocity variance showed very little dependence on the smearing parameter.
Also, in order to further evaluate the sensitivity of the simulation results to the grid resolution, the results obtained from the LSDD model with two different grids which are of sizes GR1 (480 × 96 × 192) and GR2 (320 × 64 × 128) in the x, y, and z directions, respectively, were compared. In the simulations, similar cell aspect ratios were used (∆x:∆y:∆z ≈ 2:2:1). Figure A1 also illustrates the lateral profiles of the normalized mean velocity and streamwise velocity variance obtained from the two different grid resolutions. As can be seen in this figure, there was a very small grid resolution dependence in the profiles of the mean velocity and streamwise velocity variance, especially in the far wake region at x/D = 24, where the impact of the SGS model on the turbulence level was stronger.