Turbulence-resolving simulations of wind turbine wakes

Turbulence-resolving simulations of wind turbine wakes are presented using a high--order flow solver combined with both a standard and a novel dynamic implicit spectral vanishing viscosity (iSVV and dynamic iSVV) model to account for subgrid-scale (SGS) stresses. The numerical solutions are compared against wind tunnel measurements, which include mean velocity and turbulent intensity profiles, as well as integral rotor quantities such as power and thrust coefficients. For the standard (also termed static) case the magnitude of the spectral vanishing viscosity is selected via a heuristic analysis of the wake statistics, while in the case of the dynamic model the magnitude is adjusted both in space and time at each time step. The study focuses on examining the ability of the two approaches, standard (static) and dynamic, to accurately capture the wake features, both qualitatively and quantitatively. The results suggest that the static method can become over-dissipative when the magnitude of the spectral viscosity is increased, while the dynamic approach which adjusts the magnitude of dissipation locally is shown to be more appropriate for a non-homogeneous flow such that of a wind turbine wake.


Introduction
Understanding and modelling wake-to-turbine and wake-towake interactions within wind farms have been recognised as one of the long-term challenges in wind energy research [1]. Existing wind farm wake models vary from low-fidelity empirical and semi-empirical approaches [2,3] to more sophisticated highfidelity large-eddy simulations (LES) where the turbines are parametrised using either an actuator disc [4] or an actuator line approach [5,6]. The latter, despite being the most computationally expensive method, has recently attracted the interest of the research community and has been applied to the study of many utility-scale wind farms [7e9]. This is primarily due to the increasing computational capacity of modern high-performance computer (HPC) platforms, but also as a result of the ability of LES to resolve turbine wakes to significantly finer spatial and temporal scales and thus provide deep insights into the complexities of wake interactions. It is evident that as we move towards longer and more realistic wake simulations we need to optimise the use of existing and future HPC resources in order to obtain as much information as possible from a particular problem size (e.g. expressed in terms of computational degree of freedom count). In the context of LES, we argue that this entails two factors: first the need for subgrid-scale (SGS) models that do not act in an overdissipative manner so that larger coherent structures of turbulence can be retained, and second the use of high-order schemes with spectral or "spectral-like" accuracy which can capture more flow field details with the same degree of freedom count.
Starting with the former, the development and validation of SGS models has long been an active area of research, with numerous models having been suggested and applied to turbine wakes. A general review of the various SGS models can be found in Meneveau and Katz [10] and Sagaut [11]; whereas a more wind turbine wakes specific review (including their interactions with the atmospheric boundary layer) was recently presented by Mehta et al. [12]. Sarlak et al. [13] also studied the role of SGS modelling in predicting wake statistics in wind tunnel test configurations by examining a number of explicit SGS models, including the standard Smagorinsky (SS) model [14], the dynamic Smagorinsky (DS) model [15], as well as the dynamic mixed model of Zang et al. [16]. The authors concluded that the selection of the particular SGS modelling approach is not a primary factor in determining the overall accuracy of the obtained results. Instead they argued that a good turbine parametrisation (in their particular case the use of sufficient actuator line resolution) is more important than the SGS modelling choice, and that implicit models which employ numerically dissipative upwind schemes can provide equally good results. Mehta et al. [12] also concluded that if an adequate spatial resolution is used, the impact of the SGS model choice is nullified.
Employing high-order methods with spectral or "spectral-like" accuracy has also been a common practice in LES of wind turbine wakes and wind farms. For instance, Calaf et al. [17] and Lu and Port e-Agel [18] used formulations which employ spectral schemes in the lateral directions, together with periodic boundary conditions, and energy conserving finite-difference schemes in the vertical direction. The use of periodic boundaries for spectral methods is imperative and limits the method to simple flow configurations. More recently, first Peet et al. [19] and latterly Kleusberg et al. [20] developed turbine parametrisation models within the hp/spectral element solver Nek5000 [21]. These formulations exhibit a number of advantages compared to the spectral schemes just mentioned. Not only do they provide higher accuracy in comparison to more conventional second-or fourth-order based models, they also allow complex geometries to be used in the simulations of wind farms (e.g. uneven terrain). Nevertheless, when spectral or "spectral-like" higher-order methods are used, a certain amount of numerical dissipation is still required to stabilise the numerical solution. In the case of spectral methods, a standard procedure is to truncate the smaller scales by using the so-called 3/2 rule [22]. On the other hand, when "spectral-like" accurate methods are used, this stabilisation is achieved either by employing a spectral-vanishing viscosity (SVV) operator [23,24], to selectively add dissipation to the near grid cut-off scales, or by applying a filter-based stabilisation to remove energy from the highest modes [25]. Higher-order "upwind-biased" numerical schemes can also be used [26,27] to stabilise the numerical solution. All these methods can deliver enhanced eddy viscosities and therefore can arguably be used as SGS models for LES. The former approach that uses a SVV operator for LES of turbulent flows has been demonstrated by Karamanos and Karniadakis [23]; Kirby and Karniadakis [28]; Pasquetti [29,30]; Pasquetti et al. [31]; Severac and Serre [32]; Koal et al. [33] and Lombard et al. [34]. It should be noted that the action of the SVV operator is strictly dissipative, thus the SVV-LES formulation may also be regarded and referred to as an implicit LES (iLES) technique [11], although its formulation and action is fundamentally different from that of the MILES (Monotonically Integrated LES) approach [35], for example. More recently an implicit SVV model based on the manipulation of the second derivative of high-order compact finite-difference schemes was introduced by Lamballais et al. [36] and applied in different flow configurations by Dairay et al. [37,38] and Ioannou and Laizet [39]. Dairay et al. [38] also determined the magnitude of the required SVV in the context of isotropic turbulence using Pao's equilibrium energy spectrum. Unfortunately, the proposed methodology cannot be easily extended to nonhomogeneous turbulent flows and one can only select a value based on a trial-and-error basis. However, even an optimum selected value (on the basis of better agreement with reference data) will introduce the same amount of dissipation everywhere in the computational domain without taking into account local flow features (e.g. local turbulent kinetic energy or the filtered strainrate tensor) or non-homogeneous and transitional regions. To overcome this problem, a dynamic version of the SVV was proposed by Kirby and Karniadakis [28] and successfully applied to turbulent channel flow.
Applying the SVV approach to wind turbine wake simulations will inevitably require a careful consideration of the particular flow structures. In particular, a turbine wake exhibits three distinct regions [40]: the near-wake field which extends up to approximately a diameter downstream the rotor and is dominated by the dynamics of the blade tip vortices (laminar region), the transitional (or merging zone) in which the large well-structured vortices begin to de-stabilise and eventually break up into smaller eddies (transitional region), and finally the far-wake field in which the wake flow tends e in the absence of an external forcing (e.g. an atmospheric boundary layer) e towards a self-similar wake (turbulent region). Given this inherent non-homogeneity of the flow field, it is expected that different levels of dissipation will need to be utilised. Considering all these, we pose the following questions regarding the two SVV (static and dynamic) approaches: Is the static SVV-LES appropriate for LES of high-Reynolds number wind turbine wakes? How sensitive are results to the selection of the static SVV magnitude? Can the dynamic SVV approach yield better results compared to the static approach, on the basis that dissipation will be scaled locally?
To answer these questions we make use of sixth-order compact finite-difference schemes, the implicit SVV approach of Dairay et al. [38] and a newly developed implicit dynamic SVV model which scales the amount of local dissipation with the magnitude of the strain-rate tensor, by invoking the linearity of the dissipation operator developed by Lamballais et al. [36]. As reference data against which we compare our numerical simulations we select the "blind tests" database of Krogstad and Eriksen [41] and Pierella et al. [42]. These wind tunnel measurements have been used by a number of researchers [20,43,44] and therefore represent an excellent benchmark case to compare against previous as well as future studies. To parametrise the wind turbines we make use of the actuator line technique [45] based on a turbine parametrisation which has been previously shown to better capture the key features of the near wake field (e.g. tip vortices).
This remainder of the paper begins with a presentation of the numerical discretisation and modelling techniques employed, including details on explicit SGS models, the implicit static and dynamic SVV methods and a short description of the actuator line turbine parametrisation. Subsequently, the wind tunnel tests are briefly described in section 3 and the experimental wake data are compared with the numerical model solutions in section 4. Further qualitative and quantitative comparisons and observations for the obtained wake solutions and their sensitivity to modelling parameters are also presented in section 5. A comparison between the standard Smagorinsky model and the two proposed models (static and dynamic SVV) is also given in section 6 for completeness. Finally, a summary of the main findings are presented in section 7 along with the proposed future work.

Numerical solver
To resolve the wind turbine wakes we consider the unsteady, incompressible, Navier-Stokes equations, where u i for i2f1; 2; 3g represents the velocity component in the ðx; y; zÞ directions respectively, p is the pressure field, t ij are the subgrid-scale (SGS) stresses associated with the use of an explicit SGS model, F T is the turbine actuator forcing and n, r are the kinematic viscosity and the fluid density respectively. In the implicit SVV-LES model the combined effects of the viscous term and the divergence of the subgrid-scale stresses will be represented by the dissipation operator D . As was mentioned earlier, the dissipation operator D is evaluated by taking advantage of the discretisation error of the second derivatives when high-order compact finite-differences schemes are used; details will be provided for both the standard and the dynamic cases in sections 2.2 and 2.3 below. With the adoption of the implicit SVV-LES technique, the momentum equation is modified to include the dissipation term and exclude any explicit SGS terms, It should be noted that the viscous term is also accounted for implicitly by the proposed dissipation operator as, where n s ðk; k c Þ is a wavenumber-dependent eddy-viscosity, k is the velocity-based wavenumber and k c ¼ p=D is the cut-off wavenumber associated with the mesh size D.
To solve the governing equations (1a) (or (2)) and (1b), the highorder finite-difference flow solver Incompact3d [46,47] is used. It is based on a Cartesian mesh in a half-staggered arrangement (the same mesh is used for the three velocity components u x ;u y ;u z , with a different mesh used for pressure p), sixth-order compact finitedifferences schemes [48] for the spatial discretisation, an explicit third-order Runge-Kutta method for the time advancement, and a direct spectral solver for the Poisson equation. The staggered mesh arrangement is essential in order to avoid spurious pressure oscillations on the mesh, while the compact sixth-order finite-difference schemes are essential to obtain the required "spectral-like" behaviour. To appreciate the importance of compact sixth-order accurate schemes, we briefly present in Fig. 1 the modified wave number of the first derivative and compare it to other standard central difference schemes (CDS) commonly used in energy conservative numerical solvers.
The modified wave number k 0 is obtained by conducting a Fourier analysis of the discrete first derivative of any velocity component u j , j2fx; y; zg, to obtain d is the imaginary unit and d ð/Þ denotes the Fourier transform of the function. From Fig. 1 we may observe that from all the presented schemes the sixth-order compact finite-differences scheme better matches the exact solution (k 0 ¼ k) as it remains accurate for a relative wavenumber k=k c of up to around 2=3. We refer to this scheme as a "spectral-like" accurate scheme due to its ability to capture finer spatial scales. We take care of the other discretisation errors such as the aliasing error appearing in the numerical solution by using the skew-symmetric form of the convective term as suggested by Kravchenko and Moin [49]. Finally, Incompact3d is parallelised with the aid of MPI and an efficient 2D pencil domain decomposition approach [47]. The turbine momentum source, F T , which is computed at each time step, uses a standard actuator line implementation, a description of which is briefly presented in section 2.4.

Explicit SGS modelling
To allow for a comparison between the proposed static and dynamic SVV-LES techniques with traditional explicit SGS models, we briefly present the formulation of the standard Smagorinsky model [14], an SGS model which has been very popular in the study of many turbulent flow configurations including wind turbine wakes [7,50]. Within the standard Smagorinsky model, the subgridstresses are represented using an eddy-viscosity n t as: is the strain rate tensor and The coefficient C S , also called the Smagorinsky coefficient, takes a constant value which depends on the flow configuration. For the simulations presented below we have selected a value of C S ¼ 0.16, similar to Martínez-Tossas et al. [50].

The "standard" (static) iSVV approach
Since the dissipation operator D is a non-standard one in the literature, it is instructive to discuss its implementation in greater detail. Starting with a general representation of compact finitedifferences schemes, the second derivative can be described through a 3e9 stencil formulation: where The five coefficients ða; a; b; c; dÞ can be chosen in such a way so as to ensure up to tenth-order accuracy, by satisfying exactly the following five relations Requiring only sixth-order spatial accuracy, the last two equations (7d) and (7e) do not need to be satisfied and therefore two out of the five coefficients can be chosen freely. Lele [48] showed that by choosing (a;a;b;c;d )¼(2/11, 12/11, 3/11, 0, 0) an "optimal" sixthorder accurate scheme with spectral-like behaviour can be obtained. To better appreciate the spectral behaviour of this "optimal" sixth-order scheme we may again present the modified (effective) wavenumber of the discrete second derivative and compare the modified (effective) second derivative wavenumber k 00 to the exact solution k 2 (see Fig. 2). Here we also present an expression of the modified second derivative wavenumber in terms of the five coefficients ða; a; b; c; dÞ associated with the compact scheme: Again, an excellent agreement can be observed up to approximately k=k c ¼ 2=3. The rightmost region, k=k c > 2=3 (region of small scales), is said to exhibit under-dissipative behaviour for the second derivative due to the inability of the "optimal" sixth-order scheme to effectively remove the near-grid-size scales when considered in the context of direct numerical simulations (DNS).
On the other hand, the ability to freely choose two of the coefficients in the discretisation presents a unique opportunity to add two constraints (equations) in an effort to add extra numerical dissipation and turn the discrete second derivative approximation into an SVV-like operator. The two additional constraints are chosen in such a way so that the modified wavenumber k 00 at k ¼ k m ¼ 2=3k c and k ¼ k c should match the respective wavenumber of a combined viscous term/SVV operator kernel: Inherently, by setting n 0 =n ¼ 0, the additional numerical dissipation effects (the ones that mimic the SVV operator) will be nullified and therefore a nearly exact solution may be obtained (see Fig. 2). Nevertheless, equations (9a) and (9b) together with the original equations for sixth-order accuracy (7a)-(7c) make the system complete and the five coefficients can now be expressed in terms of k The newly obtained operator is expected to exhibit enhanced levels of numerical dissipation for the higher wavenumbers due to the imposed conditions (equations (9a) and (9b)). However, to appreciate the effective spectral eddy viscosity we may assume that the modified wave number will act as any other spectral eddy viscosity model [51,52] and therefore we may argue that in the spectral space By re-arranging equation (11), we compute the implicit (or associated) spectral eddy viscosity as It can be seen in Fig. 3a that the implicit SVV closely follows the explicit SVV of Karamanos and Karniadakis [23]; which further justifies its name. The term 'implicit' refers to the fact that the resulting operator is entirely due to the controlled numerical dissipation error. In addition, when the iSVV operator is used in the context of LES (not just to stabilise the solution), we may also refer to it as the iSVV-LES model. Another important observation which can be made for the same figure is that the "optimal" sixth-order scheme attains negative values near the cutoff, a fact which again proves its under-dissipative nature. Since we previously introduced the idea of a spectral eddy-viscosity model [52], it is tempting to compare the physical-space constructed spectral viscosity n s ðkÞ with existing spectral eddy-viscosity models such as that of Chollet and Lesieur [51] which here we have scaled with the kinematic viscosity n to attain a maximum value of 1, n s ðk; k c Þ ¼ nK À3=2 where K 0 ¼ 1:114. The normalised Chollet and Lesieur [51] spectral eddy viscosity is also plotted in Fig. 3 and it is evident that the most significant difference is the low wavenumber "plateau region" of the spectral eddy viscosity which is missing from the SVV kernels and allows the spectral eddy viscosity to also affect larger flow scales. Next, looking at Fig. 3b it can be observed that the shape of the normalised spectral viscosities n s ðk; k c Þ=n 0 is largely invariant to the SVV magnitude, a property which allows us to make the following approximation for the dissipation operator: D ðu; an 0 =nÞ ¼ aD ðu; n 0 =nÞ for any an 0 =n; n 0 =n > 10: This property, will become extremely useful in the formulation of the dynamic iSVV approach to follow, as the dissipation operator will be scaled by a local quantity. Nevertheless, the magnitude of the SVV n 0 =n depends on a number of parameters including the mesh resolution, Reynolds number Re, as well as the accuracy of the underlying discretisation schemes. Karamanos and Karniadakis [23] performed SVV-LES of a smooth wall turbulent channel flow using four variations of the SVV kernel parameters, including its shape and magnitude, and compared with reference data. Unfortunately, from their simulations they were not able to identify a particular trend, although they suggested that an optimal value should exist and it should be sought through higher resolution simulations and a more meticulous comparison with the flow statistics. Later, Kirby and Karniadakis [28] obtained optimal values through comparison with reference data and higher-resolution computations. The most rigorous study to date for the estimation of the optimal value for n 0 =n can be found in Dairay et al. [38] who used a Pao-like solution for the energy spectrum [53] and were able to systematically calculate the optimal value, but unfortunately only for the case of isotropic turbulence.

The dynamic iSVV approach
In the previous paragraph, we discussed the difficulties associated with selecting the right value for the magnitude of the spectral vanishing viscosity and the fact that most studies, except when isotropic turbulence is considered, have relied on choosing their values arbitrarily. However, what is also important to notice is that by considering a uniform distribution for the magnitude of the spectral eddy viscosity within the computational domain, the SVV-LES approach cannot take into account any localised flow features, and therefore dissipation is added irrespectively. Karamanos and Karniadakis [23] recognised this defect of the static method in their concluding remarks and later Kirby and Karniadakis [28] introduced and applied a dynamic SVV model, first to the inviscid Burgers equation and subsequently, to the compressible Navier-Stokes equations. The key idea underlying the dynamic SVV approach is to scale the amount of spectral viscosity by a localised flow diagnostic. Kirby and Karniadakis [28] selected the magnitude of the local and normalised the strain-rate tensor magnitude by a global scaling quantity S ∞ ¼ maxfS g taken over the entire domain. Here we have adopted the same scaling and we propose a new implementation of the implicit dynamic SVV model by invoking the linearity properties of the dissipation operator D (equation (15)), i.e.
This approximation holds true only for large values of n 0 =n, but allows us to implement the dynamic method without changing the numerical scheme. An alternative approach would need to interact with the numerical scheme by defining a local value for n 0 =n.
Another important feature of our proposed implementation is that the ratio S =S ∞ can be calculated and monitored a priori so that additional constraints can be added, e.g. in the form of minimum or maximum values.

Actuator line model
For the turbine parametrisation, the actuator line model implementation of Deskos et al. [54] and Deskos and Piggott [55] is used. In particular, the rotor blades are parametrised into discrete blade elements and the normal and chordwise airload coefficients C n , C c are computing using tabulated lift and drag coefficients of each blade element similar to Sørensen and Shen [45] and Troldborg et al. [56]. For the presented simulations, the lift and drag coefficients were obtained from the measurements provided in the appendix of Sarlak [43]. In addition the model accounts for dynamic stall effects, using a low-Mach number modification of the Leishman and Beddoes [57] model proposed by Sheng et al. [58]; tip-loss correction effects [59], and a dynamic tower-shadow model similar to Sarlak et al. [13]. Finally, a standard Gaussian smoothing kernel approach is used to project the actuator line forces onto the fluid mesh, using a fixed smoothing parameter ε ¼ 2:2D, where D ¼ ðDxDyDzÞ 1=3 . The selection of the smoothing parameter follows the recommendation of previous studies [13,56] for uniform Cartesian meshes for which the condition 2D < ε < 4D is required.

Simulation set-up
The numerical simulation set-ups are based on the wind tunnel experiments of Krogstad and Eriksen [41] and Pierella et al. [42]. For brevity here we only present information we feel is important to maintain clarity; further details on turbine blade characteristics or flow measurements techniques can be found in Krogstad et al. [60] and Krogstad and Eriksen [41]; Pierella et al. [42]. The computational domain includes the full length of the wind tunnel which extends L ¼ 11.15 m in length, and has a width and height of W ¼ 2.7 m and H ¼ 1.8 m respectively. Krogstad and Eriksen [41] reported that their measured incident velocity profile was found to be nearly uniform with U ∞ ¼ 10 m s À1 and a turbulent intensity (streamwise velocity fluctuations) as low as I ¼ 0.3%. To reproduce such conditions in our simulations, the inlet velocity is designed from inflow planes generated by the isotropic synthetic fluctuation method of Davidson [61] using 150 modes, with lengthscales from twice the grid size to the rotors' average diameter D and a background kinetic energy k ¼ 0.01 m 2 /s 2 . Before running the turbine simulations, inflow planes are introduced and convected downstream of an empty (no wind turbines) computational domain in order to ensure that the turbulence intensity converges (after an initial dissipation of the intensity of the synthetic turbulent structures) to the required value of I ¼ 0.3%. Standard 1D convective boundary conditions are used as outflow conditions in the streamwise direction and free-slip boundary conditions are implemented in the lateral directions. All simulations presented hereafter assume an effective rotor Reynolds number of Re D ¼ 100000. For the turbine configurations we consider two cases: (1) a single turbine operating with a prescribed constant tip-speed ratio l ¼ 6 which we will refer as the "blind test 1" or BT1 hereafter, and (2) two similar turbines operating in line, with the second (rear) one being placed at a small distance (three diameters) behind the front one, which we will refer as the "blind test 2" or BT2. In the latter case the two turbines operate with prescribed rotational velocities which are different and in particular with tip speed ratios of l 1 ¼ U 1 R=U ∞ ¼ 6 and l 2 ¼ U 2 R=U ∞ ¼ 4, where U 1 and U 2 are the rotor low speed angular velocities of the front and the rear turbines, respectively. A schematic representation of the two configurations is shown in Fig. 4. The computational domain is discretised with a uniformly-spaced mesh consisting of 1281 Â 241 Â 241 nodes, and a time step of Dt ¼ 0.0005 s is used throughout.
Each blade/actuator line is discretised with 70 blade elements which ensures, together with the magnitude of Dt, stability and convergence for the actuator line/fluid solver coupling. The simulations were run for a total of 40000 time steps each and wake statistics were collected after a spin-up time of 10 s (20000 time steps). Each of the simulations was performed with 576 computational cores and ran for approximately 24 h on ARCHER, the UK supercomputing facility. It should be noted that this wall clock time could easily be reduced by increasing the number of computational cores when running the simulations.

Static SVV wake solutions
Since a rigorous approach for the calculation of the required SVV magnitude n 0 =n does not yet exist for non-homogeneous turbulent flows, we have selected the following four values: 10, 100, 1000 and 2000 through an "order-of-magnitude" approach.
The wake solutions obtained with the static implicit SVV method are compared against the reference wind tunnel data for both BT1 and BT2 in Fig. 5. There we have plotted the horizontal time-averaged normalised streamwise velocity, u x =U ∞ , and the turbulence kinetic energy, TKE=U 2 ∞ at three locations downstream the rotor for BT1, and the horizontal time-averaged normalised streamwise velocity, u x =U ∞ , and streamwise Reynolds stresses ∞ at three locations downstream the rear turbine for BT2. All four wake solutions show good overall agreement with the wind tunnel measurements, particularly for the first-order statistics (mean velocity/velocity deficit). Many wake characteristics, including the near centreline asymmetry introduced by the tower wake, the width of the wake, as well as the magnitude of the velocity deficit are all very well predicted by the model, irrespective of the value of SVV. On the other hand, the magnitude of the SVV appears to have a more pronounced effect on the turbulent kinetic energy and streamwise Reynolds stresses profiles for BT1 and BT2, respectively. Starting with BT1, the low values n 0 =n ¼ 10 and 100 are shown to better capture the near wake (x ¼ 1D) TKE and in particular the magnitude of the TKE peak on the wake's outer thin shear layer, although the discrepancy between all wake solutions and the reference data remains high. This is related primarily to the ALM turbine parametrisation, and more particularly on the choice of a uniform smoothing value for ε. It has been demonstrated that alternative formulations for the AL to mesh interpolation can mitigate these effects [62]. Nevertheless, it can be observed that an increase in n 0 =n will further smear out the two peaks. Paradoxically, an increase in the value of n 0 =n has an opposite effect on the near centreline TKE. This is clearly an effect of the linearity of the SVV operator. An increased amount of dissipation will indiscriminately annihilate smaller turbulent scales even when they experience different levels of shear and therefore create a less distributed TKE profile. In other words, the computed sharp TKE distribution is a result of "filtering" the smaller scales and allowing only a large vortex to be shed from the tower. Again, the experimental data suggests that n 0 =n ¼ 10 and 100 are much more appropriate values. Looking at the downstream profiles of the TKE and more particularly at x ¼ 3D, the difference between the four simulations appears to be more distinguishable. This is a particular region in the spatial evolution of the wake where the most energetic structures are destabilised and break down into smaller eddies. Inherently, an increase in the level of SVV dissipation may also affect the evolution of this process. The wake solution for n 0 =n ¼ 100 and 1000 yield the best match with the experimental data, which suggests that the  optimal value should lie within these two values. On the other hand by applying a larger SVV value (n 0 =n ¼ 2000) the solutions becomes over-dissipative, and lower levels of TKE are predicted. Finally, the wake solutions at x ¼ 5D (onset of far wake field) appear to converge to the same solution. This result shows that variations in the SVV magnitude become less important when the turbulent wake field becomes more isotropic. In BT2 the obtained wake statistics from the four different SVV parameters exhibit a very similar behaviour, although now the larger discrepancies appear at the profile immediately after the rear turbine (x ¼ 1D). This is due to the increased levels of TKE produced by the interaction of the upstream wake with the rear turbine, which manifests itself by first accelerating the breakup of the front turbine helical vortices, and subsequently shifting the merging/mixing zone upstream. It is worth noting that the solution for n 0 =n ¼ 2000, although quantitatively closer to the reference data, is different from the other three solutions presented in Fig. 5d. From all these observations, we may conclude that an optimal value should lie between n 0 =n ¼ 100 and n 0 =n ¼ 1000, although it cannot be established that this optimal value will out-perform all other values in terms of a rigorous quantitative metric. As we have previously mentioned, the linearity of the SVV approach will in general ignore the strength of the local flow features and will blindly add viscosity to the solution. Given the non-homogeneity of simulated turbine wake this might mean that for any selected value of the static SVV certain features of the wake, e.g. the merging zone, may be well captured whereas others may not (see Fig. 5).

Dynamic SVV wake solutions
Turning to the simulations using the dynamic SVV approach, we have chosen to run simulations only for the two higher values of SVV, n 0 =n ¼ 1000 and 2000. This decision was driven by the current formulation of the scaling factor S =S ∞ which is bound between 0 and 1, and therefore the effective magnitude of the SVV will be scaled down to S S ∞ n0 n < 10 in places, which can potentially violate condition (15). To guarantee that the local dynamic SVV magnitude does not go below 10 an extra condition has been imposed. This condition has a minimal impact on our results, and the amount of minimum viscosity can be seen as a stabilisation for the solution. Nevertheless, it is worth noting the quantitative differences between the static and dynamic SVV approaches. An immediate observation is again that the two approaches agree well as far as the mean streamwise profiles are concerned, although the dynamic method appears to better capture the wake width, particularly in profile x ¼ 3D of BT1. What is significant, however, is that with the dynamic method the discrepancies between the two SVV simulations have been significantly reduced. In particular, the TKE in the second profile (x ¼ 3D) of BT1 (see Fig. 6) does not experience the sharp drop in the merging zone peaks as we change the SVV magnitude between n 0 =n ¼ 1000 and n 0 =n ¼ 2000 as happened before, and in addition the width of the shear layer (merging zone of the wake) is better captured.

Integral rotor quantities
Lastly, from the comparison with the wind tunnel data, we consider the rotor integrated quantities such as the power and thrust coefficients. These are not expected to differ significantly for either the static or the dynamic cases, or when varying SVV parameters. These is due to the fact that the SVV operator does not affect the large energetic scales. Nevertheless, the power and thrust coefficient are defined via, (17a) where Q is the generated shaft torque, U the prescribed rotational velocity of the rotor and A ¼ pR 2 the rotor's swept area. We compute the time-averaged power and thrust coefficients after a spin-up time of 20 s which corresponds to at least 40 revolutions and the computed values are compared against the respective experimental values of Krogstad and Eriksen [41] and Pierella et al. [42]. The final relative error estimates for all the simulations are shown in Table 1.
It is shown that all simulations exhibit similar levels for the relative errors for both the power and thrust coefficients. More specifically, in all simulations including both BT1 and BT2, the estimated power coefficients are over-predicted by approximately 3% for BT1 and 5% for the BT2 front turbine and 6.5% for the rear one. On the other hand, C T is under-predicted for BT1 and the BT2 front turbine (8e10%), while slightly under-predicted for the rear BT2 one. Nevertheless, it is worth noting that for BT2 a bigger scatter is observed for the rear turbine quantities. This is not a surprising effect as the second turbine is directly affected by the strength and turbulence structure of the upstream wake. We observe that the dynamic SVV for n 0 =n ¼ 1000 provides the smallest error for both C P and C T for the rear turbine, while by increasing the dynamic SVV magnitude to n 0 =n ¼ 2000, the error is increased by less that 1% for both. Finally the values for BT1 and the front turbine of BT2, although they only experience inflow isotropic turbulence, also exhibit a better match with experimental values and less variability when the dynamic SVV method is used. This might be attributed to the interaction of the near wake field with the tower wake which the rotor experiences as dynamic inflow.

One-dimensional spectra and wake visualisations
One of the main advantages of LES is that it can capture the large coherent structures of turbulence with high accuracy. As we mentioned in the introduction, wind turbine wakes are wellknown for their near-wake field helical vortices, as well as the particular mechanism which leads to them breaking up into smaller eddies. Stability analysis of the near-wake helical vortices have been undertaken both in the context of LES [44,63,64] and analytical models [65e67]. In all these studies the main mechanism that leads to an effective breakup of the helical structures has been identified to be the merging and pairing of the neighbouring vortices caused by small perturbations present in the ambient flow.
In an effort to highlight the ability of turbulence-resolving SVV simulations to capture the aforementioned mechanism as well as to demonstrate the advantages of the dynamic SVV, we present onedimensional velocity spectra and wake visualisations for qualitative comparisons.

One-dimensional spectra
A better appreciation of the dissipation of the two methods (static and dynamic) as well as their respective selected parameters, can be sought through consideration of one-dimensional energy spectra. To obtain a meaningful estimate for the energy spectrum of the wake, we have deployed a total number of 30 "probes" (15 along the centreline and 15 along the periphery of the  wake, placed within half a meter distance each) downstream the turbines to sample the instantaneous velocities. Applying the frozen turbulence theory of Taylor [68,69], we estimate the power spectrum density (PSD) for each component of the fluctuating velocity (u 0 i ¼ u i À u i ) which are presented together with the characteristic inertial range k À5=3 slope of Kolmogorov [70]; referred to as K41 hereafter in Fig. 7a and b for BT1 and BT2, respectively.
The dimensionless power density spectra appear to follow similar trends irrespective of the selection of either the SVV magnitude for both BT1 and BT2. In particular, in the lower frequency regime the different power spectral densities remain constant (large energetic scales) while subsequently all wake solutions follow the inertial sub-range À5=3 slope of the Kolmogorov spectrum [70]. Perhaps the only region in which important differences can be observed is at the right-most regime. This is not surprising as the implicit SVV model injects enhanced levels of dissipation according to the magnitude of SVV at the smaller scales and therefore it should affect the "tail" of the energy spectrum. Taking into account that the shape of the spectrum's tail is directly connected to the effective filter used in LES [11], the differences that are observed in all solutions confirm the ability of the implicit SVV-LES method to act as an explicit SGS model. Moreover, the PSD obtained from the dynamic SVV solutions can better capture the energy distribution along the relevant scales, as it lies within the two curves n 0 = n ¼ 10 and n 0 =n ¼ 1000 in the right-most region. The exact effects of filtering may be obtained by computing the transfer function, where E SVV ðkÞ is the wavenumber-dependent energy spectrum of the SVV simulations and E DNS ðkÞ is the direct numerical simulation (DNS) respective one. Unfortunately, the high Reynolds number assumed for the simulations makes the computational cost of the respective DNS simulations prohibitively large.

Wake visualisation
A common way to identify the coherent structures of a fluid flow is by using the so-called Q-criterion of Hunt et al. [71]; which describes the vortical structures via the positive second invariant of the velocity gradient. A similar technique was introduced by Jeong and Hussain [72] with the so-called l 2 -criterion (not to be confused with the definition of the tip speed ratio used in the present work's nomenclature). To generate results comparable to previous studies [13,44] observe the qualitative differences between the low dissipation and the high dissipation solutions. To appreciate the quality of the plotted coherent structures, it is worth re-iterating the mechanism of the near-wake destabilisation and evolution to the far-wake field. Instabilities start to appear at the end of the near-wake field and manifest themselves via intermediate-wave modes, which are multiples of the dominant tip vortex frequencies, and evolving along the tip vortices. Widnall [67] showed that such instabilities are due to the mutual-inductance of the tip vortices and that they have many common features with the vortex pair instability. Vortex pairing together with a subsequent vortex merging mechanism as was later demonstrated by Sarmast et al. [64] is responsible for the destabilisation of the near wake and its transition to the far wake field. We should emphasize at this point the importance of using a synthetic inlet method in triggering these modes and invoking the so-called "pairing-and-merging" mechanism. Numerical experiments with grid-scale randomly generated inflow noise or uniform inflow profiles were found to be unable to yield the desired results.
Looking at the obtained vorticity iso-contours plotted in Figs. 8 and 9, we may observe that the above mentioned mechanism is visible only in the instantaneous coherent structures obtained with the static SVV approach for n 0 =n ¼ 100 and 1000 as well as the dynamic SVV approach for 1000 and 2000. All other solutions do not capture these features satisfactory. For instance, the solutions obtained by using small amounts of dissipation (e.g. n 0 =n ¼ 10) are found to exhibit short-wave instabilities, which most likely are a result of grid-scale spurious numerical oscillations. On the other hand, in SVV solutions with large amounts of numerical dissipation, the effective spectral viscosity would interact with the mutualinductance modes/scales and therefore it will have an impact on the transition mechanism. This is confirmed by the wake visual-   disappeared from the iso-contours. Apart from the transitional mechanism, the far wake field is also affected by increasing the magnitude of the SVV, as only the larger coherent structures remain in the solution. Thankfully, due to increasing levels of isotropy in the far wake field, an increase in the SVV magnitude will not affect any vital mechanism but instead it can be seen solely as an increase in the effective filtering as noted by Dairay et al. [38]. A "golden ratio" type of solution can be obtained by the dynamic SVV approach (e.g. DSVV with n 0 =n ¼ 1000), as the local scaling of the dissipation acts more appropriately for both the near/transitional and far wake fields.

Comparison with the standard Smagorinsky model
Having discussed the quantitative and qualitative differences between the static and the dynamic SVV methods as well, as their sensitivity to the modelling parameter n 0 =n, it is instructive to make a comparison with the solutions obtained using the standard Smagorinsky model. For brevity, only results for the BT1 case are shown. In Fig. 10a and b, three wake solutions are plotted together with the corresponding experimental data. These are the static SVV model with n 0 =n ¼ 1000, the dynamic SVV model with n 0 =n ¼ 1000 and the standard Smagorinsky model with C S ¼ 0.16. All other parametrisations including the turbine forcing and inlet turbulence were kept the same. The obtained solutions appear to be very similar for the mean streamwise velocities for all the profiles downstream of the rotor. However, the effects of the different SGS modelling techniques appear to be more pronounced in the second-order statistics (TKE). More precisely, it may be observed in the TKE profiles that the two peaks corresponding to the mergingzone wake solution are under-predicted by the Smagorinsky model compared to the dynamic SVV model. This can be attributed to the overly dissipative action of the Smagorinsky model for the larger and intermediate turbulent scales (no wavenumber-dependence for the eddy-viscosity). On the other hand, the static SVV model appears to be less accurate than the Smagorinsky model in the vicinity of the near-wake field, and in particular at the wake interface with the ambient fluid. Evidently, the static SVV model dissipates the smaller-scales of the near wake field irrespective of the shear-rate tensor magnitude which results in over-damping of the outer thin shear layer. On the other hand, the Smagorinsky model which takes into account the local shear-rate tensor does not act as dissipatively. Ultimately, at x ¼ 5D, the three solutions appear to converge although differences can still be observed and both the static and dynamic SVV methods are a better match with the experimental data. All in all, the dynamic SVV method appears to better capture the dynamics of the non-homogeneous, anisotropic flow field as it combines a local adjustment of dissipation (via the shear-rate tensor magnitude) with a scale-dependent (wavenumber-dependent) eddy-viscosity model.

Discussion and conclusions
We have investigated the ability of static and dynamic spectral vanishing viscosity (SVV) operators to act as a SGS model and predict the wake statistics behind a single wind turbine as well as two turbines operating in-line. This was motivated by previous studies in the field of wind turbine wakes [12,13] which have pointed out that implicit LES formulations which rely on numerical dissipation for the SGS modelling can predict the wake equally well. In the present formulation we consider a higher-order compact finite-differences flow solver, the implicit SVV approach of Dairay et al. [38]; and a newly developed implicit dynamic SVV approach to conduct LES of wind turbine wakes. The accuracy of the implicit SVV approach (static and dynamic) relies on the selection of the SVV magnitude and therefore a parametric analysis was required to understand its behaviour. To appreciate the sensitivity of the two approaches to the input parameter as well as determine an optimum value we perform comparisons against reference wind tunnel wake data. From this analysis we found that for the static SVV model an optimal value exists which provides a better match with the reference data, and this was found to lie between n 0 =n ¼ 100 and 1000. This optimum value was confirmed in both BT1 and BT2 test cases. On the other hand, the dynamic SVV showed a better match to reference data than static SVV, while a considerable variation of the n 0 =n parameter (from 1000 to 2000) did not impact their quality. A brief comparison between the standard Smagorinsky model and the two methodologies proposed in this work showed that the dynamic SVV model provides solutions which better match the experimental data. Moreover, an advantage of the implicit SVV models, as compared to other standard or more sophisticated SGS models, is that they require no extra computational cost associated with the calculation of the implicit SVV term as their action is already incorporated in the numerical schemes. A minor increase in the computational cost is required by the dynamic SVV model due to the calculation of the global shear-rate tensor magnitudes (the local ones are already computed for the advective terms of equation (2).) More specifically, we followed a heuristic approach in order to conduct a parametric analysis for the magnitude of n 0 =n which is a priori unknown to us. The values of n 0 =n were chosen with an "order-of-magnitude" based approach and the analysis showed that the optimal should lie within 100 and 1000. The linearity of the static SVV was shown to have two main drawbacks. First, as the same magnitude of dissipation was added to all the grid cut-off scales, irrespective of the local flow features (e.g. strong shear rate), even when using a nearly optimum value the produced solution created an ambiguous picture for the accuracy of the method. For instance, in the wake prediction for BT1 the strong near-wake shear layer was smeared out reducing the match with the experimental data while the downstream TKE predictions (near the merging zone region) agreed with them rather well. This means that while the value of 1000 is good for one region of the flow its performance is very poor for another. This effect was not observed in the study of Dairay et al. [38] and it is believed to be a feature of non-homogeneity as well as the interaction of the linear operator with the transitional part of the wake. The second drawback of the static method is related to the time-step stability of the simulations, particularly because an explicit third-order Runge-Kutta method was used here. Inherently, in our simulations the viscous stability limit of Lele [48] constraints the time step in order to maintain stability in the numerical solution. When the magnitude of SVV exceeds a value of 1000 the initially assumed Dt ¼ 0.0005 s did not satisfy the stability condition any more and therefore the time step had to be reduced by half. Surprisingly, the same condition was not required for the dynamic SVV approach, a finding which needs to be further examined in future studies. Nevertheless, apart from better time-step stability properties, the dynamic SVV addresses the problem of an indiscriminant dissipative action of the static SVV by scaling the local SVV magnitude with the magnitude of the local shear stress tensor. This property was shown to provide better results and rely less on the selection of the n 0 =n parameter. In fact the dynamic SVV approach has more similarities with the standard Smagorinsky model [14]. Its dependence on the local flow features has been shown to provide better results both in a qualitative and quantitative sense.
In summary, the present study has confirmed previous findings by showing that wind turbine wakes can be predicted with high accuracy if well-resolved simulations are considered, and that an implicit LES approach which depends on numerical dissipation can serve as an SGS model equally well. In the framework of SVV-LES or iSVV-LES models, the static method provides good results when an optimal magnitude value is used. However, a dynamic SVV approach would be more appropriate particularly when nonhomogeneous turbulent flows are considered. Future studies will focus on improving the representation of the local scaling parameter so that the dynamic SVV method can be easily applied to homogeneous flows as well. We should also emphasize that the term "dynamic" can potentially create confusion as to the nature of the method. The term dynamic, is usually used to point out that the method behaves in an autonomous manner by adjusting all relevant parameters in order to achieve a goal (e.g. minimise energy dissipation). Here the term dynamic is used to signify that the value of SVV is not homogeneous but varies in space and in time with the magnitude of the shear-rate tensor. A truly dynamic method should be able to use flow diagnostics to determine the SVV magnitude, similar to the dynamic model of Germano et al. [15].
Future work will focus on extending the SVV methodologies to more realistic wind energy flow configurations, including the ability to capture the turbulent structures of the atmospheric boundary layer (ABL). In such flow conditions, the selection of the SVV magnitude n 0 =n and its adjustment near the wall, are expected to have a significant impact on the final solutions.