A Numerical Study on the Development of Self-Similarity in a Wind Turbine Wake Using an Improved PseudoSpectral Large-Eddy Simulation Solver

Large-eddy simulation (LES) is performed to investigate self-similarity in a wind turbine wake flow. The turbine is represented using an actuator line model in a pseudo-spectral method-based solver. A new hybrid approach of smoothed pseudo-spectral method and finite-difference method (sPSMFDM) is proposed to alleviate the Gibbs phenomenon caused by the jump of velocity and pressure around the turbine. The LES is validated with the mean velocity and turbulence statistics obtained from wind-tunnel measurement reported in the literature. Through an appropriate choice of characteristic scales of velocity and length, self-similarity is elucidated in the normalized mean velocity and Reynolds stress profiles at various distances. The development of self-similarity is categorized into three stages based on the variation in the characteristic scales and the spanwise distribution of normalized velocity deficit. The mechanisms responsible for the transition of self-similarity stages are analyzed in detail. The findings of the flow physics obtained in this study will be useful for the modeling and fast prediction of wind turbine wake flows.


Introduction
Wind power is one of the fastest-growing sources of electricity supply, with advantages in environmental benefits and the abundance in resource [1].When designing the layout of wind farms, on one hand, it is desirable to have the distance between turbines along the wind direction as short as possible for less land use and lower cost of auxiliary facilities.On the other hand, insufficient distance between turbines leads to less wind energy extraction.The characteristics of the mean velocity in the wake of a turbine is a key factor in determining the lower bound of turbine array spacing.
The mean flow of a turbine wake is influenced by many variables of the ambient flow and the turbine structure.If the turbine rotor is simplified as a static penetrable disc perpendicular to the mean flow, classical boundary-layer theories could provide analytic predictions for the development of the free wake in laminar or turbulent flows [2,3].The analytic predictions usually start the trial solution from the knowledge in the self-similarity phenomenon, which emphasizes the similarities of the spanwise distribution at various downstream distances.The theoretical solution shows that the inflow velocity, the drag coefficient, and the Reynolds stress are the dominant factors affecting the mean flow in the axisymmetric wake.The drag coefficient of the turbine equals to the thrust coefficient, which can be predicted using the blade element momentum theory (BLMT) [4,5].Based on the BLMT, the drag coefficient is determined by the blade configuration and tip-speed-ratio.The effect of the ambient turbulence on the wake flow, which accelerates the recovery of velocity [6,7], is reflected in the Reynolds stress term in the boundary-layer equation for the wake flow.The Reynolds stress is often estimated using the mixing-length model [2,8].
In reality, the wake flow of a wind turbine is more complex than the simplified axisymmetric disc wake, owing to the rotation of the blades and the blockage of the nacelle.The rotation of the blades generates helical tip vortices and root vortices.The nacelle reduces the velocity at the center of the rotating plane and generates a secondary wake flow inside the turbine wake.In addition, there exists the phenomenon of wake meandering, which is a large-scale lateral motion of the central low velocity region under the perturbations of ambient turbulent flow.The above various features in the turbine wake may interact with each other.The interaction between tip vortices and disc wake flow has been investigated theoretically [9] and experimentally [10], where the mutual inductance between adjacent spirals is found to be the dominant mechanism for the merging and grouping of tip vortices.For some blade geometries, strong root vortices may form at a relatively large radial distance from the centerline and further interact with the tip vortices [11].The wake meandering has been studied in laboratory experiments using particle image velocimetry (PIV) [12].In numerical simulations [13], it is found that the amplitude of wake meandering is proportional to the thrust coefficient and inversely proportional to the frequency of the hub vortices.The mean velocity shear in atmosphere also makes the inflow condition different from that in the classical wake theory.Even so, the high turbulence intensity produced by the large mean velocity shear is found to have stronger effect than the high mean velocity shear itself [7].
Large-eddy simulation (LES) is a powerful and increasingly popular tool in the study of turbine wake flows.The complexity of the wake structure makes it challenging for theoretical and experimental investigation to gather the whole-field data for a wide range of parameters, while setting up cases for phenomena of interests is more flexible in LES.The actuator line method (ALM) and its predecessor actuator disk model (ADM) have been extensively applied in LES to consider the effects of turbine forces on the wake flow with a computational cost lower than the geometry-resolved LES (LES-GR) [14][15][16][17][18][19][20][21][22][23][24][25][26].While the combination of LES and ALM/ADM (LES-AL and LES-AD) has been widely employed for the estimation of the mean flow and the corresponding wind energy in wind farms [15,17,18,21,[23][24][25], less LES studies have been performed to study the mechanisms of the interactions among different flow structures in the turbine wake.Bastankhah and Porté-Agel [27] derived an analytical wake model that assumed a Gaussian distribution for the velocity deficit in the wake while considering the conservations of mass and momentum.They compared the results of the analytical model with wind-tunnel measurements and LES.The analytical model predicted the velocity deficit and power extraction more accurately than previous models that considered a top-hat shape for the velocity deficit [28,29].Xie and Archer [30] used LES-AL to study the self-similarity of the mean velocity and the added turbulence intensity in a single turbine wake flow.They found that the self-similarity of velocity deficit exists in the horizontal and upper vertical planes, while it breaks down as it gets close to the ground in the lower vertical plane.They proposed to use two different wake widths for horizontal and vertical planes, respectively.Their modified wake model provides better predictions for the mean velocity deficit and turbulence intensity in far wake than previous models [27][28][29].Kang et al. [19] compared the results of LES-GR and LES-AL, and emphasized the necessity of nacelle modeling in simulations to capture the radial expansion of hub vortex meandering.Foti et al. [31] found that the wavelength of the wake meandering profile becomes larger after the expanding hub vortex intercepts the tip shear layer.These above studies provided quantitative analyses focusing on nacelle model, tip vortices, wake meandering, or specific aspects of the self-similarity.In the present study, the interaction of multiple flow structures in the wake is investigated.
The present work aims to use LES-AL to study the development of self-similarity in the wake of a turbine, including the formation of self-similarity and the influences from dominant flow structures in the wake.The rest of this paper is organized as follows.In Section 2, we provide an overview of the computational framework, which embeds the actuator line model into a spectral solver of LES.A new spectral scheme for calculating derivatives is proposed to alleviate the Gibbs phenomenon observed in the coupling of ALM and spectral solver.In Section 3, we set up a simulation for a miniature turbine in a wind-tunnel test reported in literature.In Section 4.1, we compare the mean velocity, turbulence intensity, and Reynolds stress at two different downstream locations with the wind-tunnel measurements.In Section 4.2, we study the self-similarity of the mean velocity by normalizing the flow field with characteristic scales of velocity and length.In Section 4.3, we present the self-similarity in the Reynolds stress.In Section 4.4, we explain the formation of self-similarity by a budget analysis of the mean momentum equation.In Section 4.5, we analyze the influences of tip vortices and wake meandering, which are assumed to cause the change in the growth rate of the self-similarity of the mean velocity.The conclusions of this study are provided in Section 5.

Large-Eddy Simulation of Turbulent Wind Field
In this study, a neutrally stratified atmospheric boundary-layer flow is considered for the wind field around a wind turbine.In LES, the fluid flow is described by the filtered Navier-Stokes equations for incompressible flows as follows, The coordinates are denoted as x i (i = 1, 2, 3) = (x, y, z), where x and y are respectively the streamwise and spanwise coordinates and z is the vertical coordinate, with z = 0 corresponding to the ground.The velocity components in x-, y-, z-directions are denoted as u i (i = 1, 2, 3) = (u, v, w).In Equations ( 1) and ( 2), ( • • •) indicates a spatial filtering at the grid scale ∆; ρ a is the density of air; τ ij = u i u j − ũi ũj is the subgrid-scale (SGS) stress tensor, and τ d ij is its trace-free part, which is calculated using the dynamic Smagorinsky model [32]; p * = p + τ kk /3 − p ∞ is the filtered modified pressure.An imposed pressure gradient ∂p ∞ /∂x is applied in the turbulent inflow precursor simulation and wind turbine simulation to drive the flows to statistical steady states.The turbine force f tb takes effect around the turbine grid points.
The fractional step projection method [33] is used in the simulation of the momentum and continuity equations.The momentum equations are rewritten as where the asterisk in the pressure term is omitted and F u is a collection of the right-hand side terms in Equation (1) except for the dynamic pressure gradient term, Equation ( 3) is split into two steps with the introduction of the intermediate velocity u * [33], The pressure is solved based on u * , 1 The velocity components are updated from u * with the correction of pressure gradient (Equation ( 6)).
For temporal discretization, a second-order Adams-Bashforth (AB2) scheme is employed for the convective terms and residual stress terms in Equation ( 4).For spatial discretization, in the horizontal directions x and y, a Fourier-series-based pseudo-spectral method is used on a uniform mesh.In the vertical direction z, a second-order finite-difference method is used on a staggered grid.Similar schemes for the spatial discretization are used in the atmospheric boundary-layer studies [15][16][17]26,34,35].A free-slip boundary condition is applied on the top boundary, and a wall model is applied on the bottom boundary as Here, τ is the shear stress, d 1 is the vertical size of the first off-bottom grid, z 0 is the characteristic roughness, at which the mean flow velocity is zero based on the log-law of wall, and κ ≈ 0.4 is the von Kármán constant.

Turbine Model
The actuator line model is used to calculate the turbine force exerting on the air flow.The blades are discretized into line segments from the blade roots to the tips.The aerodynamic lift force F L and drag force F D on each line element are evaluated based on the local relative flow velocity U rel , empirical drag coefficient C D and lift coefficient C L , chord length c, and its length ∆r in the radial direction of the blades: The empirical coefficients C L and C D are interpolated from their database based on the angle of attack α = β − γ, where γ is the twist angle of the airfoil section at the location of the line element, and β is the angle between the relative velocity U rel and the rotation plane, as shown in Figure 1.
The turbine is discretized on a Lagrangian mesh whose coordinate vector X is time-dependent, while the fluid is discretized on a Eulerian mesh whose coordinate vector x is time-invariant.The flow velocity u(X) at the Lagrangian coordinate X of a blade line element is evaluated from the velocity field u(x) on the background fluid grid.The aerodynamic force f (X) of a line element is the vector sum of F L and F D in Equations ( 9) and (10) and then distributed to the fluid grid as f (x).We use the discrete Dirac function to project velocity from the background fluid to the blade grid points.
where the displacement vector x − X is normalized by the grid size ∆x i in the discretized domain, and δ(r) is a smoothed version of the discrete Dirac function [18,36], To project the force from the blade location to the fluid grid, we use the smoothed Dirac function when the grid is relatively coarse, or the Gaussian kernel when the grid is fine enough to resolve the kernel width ≈ 0.17c [22,37], Here, N b is the number of blades and x − re i is the distance between the fluid grid point location x and the i-th blade element location re i , of which the direction vector is e i .The influential ranges of summations for u(X) and f (x) are limited as the delta function and the Gaussian function vanish at large distances.
The nacelle and tower are discretized by triangular surface elements and their forces on the air flow are estimated following the idea of the immersed boundary method [38], Here, the force f (X) recovers the flow velocity ũ(X) at the surface element to the desired velocity u d (X) after a timestep ∆t.The coefficient C se for all the surface elements assures the drag coefficient C D for the bulk solid to be a preset value, where U ∞ is the undisturbed inflow velocity, A is the cross-sectional area of the structure (nacelle and tower), and e is the unit vector in the inflow direction.The upstream fluid velocity is decomposed into a parallel component u and a perpendicular component u ⊥ with respect to the rotation axis.The rotation speed is the product of the angular velocity Ω and the distance r to rotation center.U rel is the relative velocity in coordinate system of the airfoil section.

Treatment of Gibbs Phenomenon in Pseudo-Spectral Solver with Turbine Model
It is known that in spectral methods, the spatial discontinuity of a variable leads to unphysical oscillations called the Gibbs phenomenon.In the simulation of turbine wake flows, the existence of turbine force may lead to spatial jump for variables of interest, such as the pressure and velocity.To alleviate the Gibbs phenomenon due to the discontinuity around the turbine, modifications are needed for the spectral method.
The turbine force f tb affects the Navier-Stokes equations as follows.There exists a jump in f tb at the turbine location, which leads to a jump in F u through Equation ( 4) and a jump in u * through Equation (5).Then, the pressure p experiences a jump through Equation (7) and finally influences the updated velocity u n+1 through Equation (6).The spectral method employed in the calculation of spatial derivatives, such as ∇p and ∇ • u * , may produce unphysical oscillations if the jump is not well treated.
To alleviate the Gibbs phenomenon in the spectral solver, we propose a four-step scheme for the calculation of derivatives: identification, smoothing, derivation, and correction (ISDC).Previously, Li et al. [39] found that the numerical error from the discontinuity can be prevented from propagation to other regions of the computational domain if the discontinuity is smoothed before transferring the values to the right-hand side terms of the Poisson equation.They tested the effect of smoothing on the problem of flow past cube, where the spectral method was combined with the immersed boundary method (IBM).We extend their idea of smoothing to the ISDC scheme for the problem of flow past wind turbine, where the spectral method and ALM are combined.Among the four steps of ISDC, the identification step is to set flags for the location of discontinuity.In our code, the rotation of blades requires the identification at each time step.The flags are used to indicate the beginning and end of each segment of the discontinuous region.The smoothing step is to remove the discontinuity, which is described with an example later.The derivation step is to employ the classical spectral method for the computation of derivatives.Because the input of spectral derivation has been smoothed, the Gibbs phenomenon is alleviated.The correction step is to employ the finite-difference method for derivatives in the discontinuous region, which ensures the accuracy there.The new scheme is applied in all gradient calculations except for the pressure gradient in the Poisson equation.For the latter, a high-order smoothing is applied after the right-hand side term is transformed from the physical domain to the spectral domain [40].
The new approach for derivative calculation is denoted as "sPSMFDM", which means a combination of the smoothed pseudo-spectral method and the finite-difference method.The procedure of this hybrid approach is illustrated in Figure 2 with a test case of the one-dimensional inviscid Burgers equation.In the region where discontinuity occurs, the variable is interpolated with a cubic polynomial that has an effect of smoothing.Two points on each side of the discontinuous region are used to determine the coefficients in the cubic polynomial.After the interpolation and smoothing, the pressure field has less discontinuity and can be safely transferred to the classical pseudo-spectral method function.The output of the pseudo-spectral method (PSM) function should be highly accurate outside the discontinuous region.For derivatives inside the discontinuous region, the finite-difference method (FDM) is applied as a correction to the output from a smoothed input.As this FDM correction is needed only for a very small portion of grid points, the increment in computational cost is negligible.Our tests indicate that in the presence of discontinuity, the sPSMFDM provides derivative estimation with higher accuracy than classical PSM.A numerical test for sPSMFDM is provided in Appendix A.

Numerical Experiment
A numerical experiment is performed in the above-described LES-AL framework.The LES-AL code is modified from the LES framework described in [18,21,41,42].The ALM part of the code is modified from the corresponding module in VFS-Wind, an open-source CFD code described in [20,43].Similar LES-AL codes are described in [16,17,26].
A miniature three-blade horizontal axis wind turbine is represented by actuator lines and placed in a numerical wind tunnel.The physical setup corresponds to the laboratory study performed by Chamorro and Porté-Agel [44], which was also adopted in other numerical studies [16,18,20,45] for method validation.The computational domain, as shown in Figure 3, has a size of L x × L y × L z = 4.32 m × 0.72 m × 0.46 m, where the height L z is chosen based on the boundary-layer thickness reported in the experiment.The boundary condition for the streamwise velocity at the top is free-slip, ∂u/∂z = 0, while the bottom boundary is non-slip, u = 0.A velocity relaxation zone [16,46] is used in the simulation to enforce an inlet boundary condition at the upstream.A precursor simulation, which develops to a statistically steady state without the turbine inside, provides a turbulent inflow of friction velocity u * = 0.102 m/s with surface roughness z 0 = 0.03 mm.A three-blade miniature turbine, which has a blade diameter of D = 0.15 m and hub height of z hub = 0.125 m, is deployed at x = 0.75 m, with x = 0 corresponding to the inlet boundary.The actuator line parameterization for the turbine is set based on the prototype rotor GWS/EP-6030×3 [44].The airfoil type along the entire blade is simplified as a circular arc airfoil.The lift and drag coefficients of the cambered airfoil at low Reynolds number are measured in the experiment by Sunada et al. [47] and provided in the simulation by Stevens et al. [26], as shown in Figure 4a.For angles of attack that are not listed in the experiment report, the aerodynamic coefficients are evaluated by theoretical formulas [48] where C f is the drag coefficient at zero attack angle, and the threshold of small and high attack angle α is adopted as the stall angle of the airfoil type.The chord length and twist angle at a normalized distance r/R to the rotation center were obtained in the numerical study of Wu et al. [16], which are plotted in Figure 4b.The actuator line elements start from the radial distance r = 0.01 m, while in the inner region r < 0.0075 m the nacelle model takes effect following Equation (15).Below the nacelle, the tower is assumed to be a vertical cylinder of diameter D tw = 0.005 m and length L tw = 0.118 m.The tower drag is considered in the same way as the nacelle drag.The drag coefficients of the nacelle and tower are 1.5 and 0.8 respectively, which are set based on the empirical ranges and our tests to acquire a best fit of the mean velocity profile in the near wake at x/D = 2 compared to the experiment measurement.A constant rotation speed Ω of the turbine rotor is determined based on a tip-speed-ratio of λ = ΩR/U hub = 4.0 and the hub mean flow velocity U hub = 2.02 m/s.The computational domain is discretized on a structural Cartesian mesh of N x × N y × N z = 384 × 96 × 65 cells.The grid number is refined from a previous simulation by Porté-Agel et al. [17] in the same physical domain, which adopted a similar spatial scheme with ours in a coarser mesh of N x × N y × N z = 192 × 32 × 42 and achieved plausible accuracies in both the mean velocity and turbulent intensity profiles.The first off-wall grid point is located at z + = z 1 u * /ν = 30.The shear stress at the bottom are obtained following Equation (8).The grid is stretched in the vertical direction so that it is denser in the rotating region.The blade diameter is covered by about 24 grid points in the vertical direction, which is refiner than 20 grid points per diameter as recommended by Yang et al. [20] for the simulation of the same experiment.Each blade is discretized into 40 segments with a resolution refiner than the fluid mesh.The nacelle and tower are represented by triangular surfaces.The grid of the fluid domain and the turbine structure are shown in Figure 5.The simulation runs with a constant time step ∆t = T/120, where T is the period of turbine rotation.After the bulk flow advects in the streamwise direction for one cycle of the domain, the simulation continues running for 235T to perform statistical analyses for the turbine wake flow.The proposed method sPSMFDM is applied in the simulation to alleviate the Gibbs phenomenon.As pointed out by Martínez-Tossas et al. [37], there is an optimal smoothing length scale for ALM, which is around 14%-25% of the chord length of the blade.This optimal smoothing length scale is smaller than our averaged grid size.To make our smoothing width closer to the optimal value, we used the three-point discrete Dirac function in Equation ( 12) which results in the Gibbs phenomenon if a standard PSM is used for the calculation of the spatial derivative.The numerical oscillations can be observed in the distribution of the spatial derivatives, such as the streamwise advection of the turbulent kinetic energy (TKE, k = u i u i /2) shown in Figure 6a.After applying the proposed special numerical treatments in Section 2.3, the numerical oscillations are alleviated in the wake region of the domain, as shown in Figure 6b.It should be noted that there are three types of spectral operators in the code, namely the spatial derivative operator, the dealiasing operator, and the spectral solver for pressure.The special treatments introduced in Section 2.3 take care of the first and the second spectral operators, while the spectral solver for pressure is unchanged.Therefore, instead of totally eliminating the numerical oscillations, the special treatments in this study alleviate the numerical oscillations and restrain the remained oscillations to a smaller region in the very near wake.The treatment for the spectral pressure solver is to be studied in future.

Validation of the LES Coupled with the Actuator Line Model and Nacelle Model
The time-averaged profiles of velocity, turbulence intensity, and Reynolds stress at 2D and 5D downstream of the turbine are shown in Figure 7.In Figure 7a,d, the mean velocity profiles of our simulation agree well with the wind-tunnel measurement [44].The influence of turbine forces is reflected by the mean velocity deficit as compared to the undisturbed upstream flow.The initial deficit in the vicinity of the turbine is caused by the energy extraction of the blades and the drag of subsidiary structures such as the nacelle and tower.As the distance to the turbine increases, the velocity deficit decays.In Figure 7b,e, the turbulence intensity profiles of our simulation show similar trend as the measurement data.The magnitude of turbulence intensity in simulation is smaller than that of wind-tunnel measurement, which is not surprising as the large-eddy simulation does not resolve all scales of the turbulence [25].The peak of turbulence intensity at the rotor top contains the contributions of tip vortex and large vertical gradient of the mean velocity, while the turbulence intensity at the rotor bottom is less or equal to that of the inflow due to negative or small velocity gradient [49].In Figure 7c,f, the Reynolds stress profiles of our simulation also show similar trend as the measurement data.The simulation profiles of turbulent components are closer to the measurement at x/D = 5 than at x/D = 2.In summary, the simulation predicts the mean velocity accurately and reproduces the turbulence properties reasonably well.In the following sections, more flow characteristics are presented from the simulation data.

Self-Similarity in Mean Velocity
Self-similarity is widely observed in free shear flows such as mixing layer and plane jet.Because the rotation region of turbine blades is a circular disk, the wake of turbine may conform to the theory of axisymmetric wake.In the vertical cross-section, the planar boundary layer where the turbine sits introduces asymmetry to the wake flow.On the other hand, the wake flow in the horizontal cross-section remains symmetry along the centerline, as shown in Figure 8a.
In classical boundary-layer theories, the mean velocity in far wake is assumed to satisfy the trial function [2] U(x, y) where U N (x) is the characteristic velocity difference.The original meaning for the subscript "N" might be "negative".In this study, "N" indicates that the variable is a characteristic scale used for "normalization".∆ is the characteristic width of the spanwise distribution of velocity deficit, and η is the non-dimensional spanwise distance.The trial function separates the velocity deficit U(x, y) − U ∞ to contributions from streamwise distance x and normalized spanwise distance η, respectively.For a turbulent axisymmetric wake, the spanwise distribution of velocity deficit is close to the Gaussian distribution, as shown in Figure 8b.The maximum deficit |U N (x)| at each streamwise distance x is chosen as the characteristic velocity difference.The half width where the velocity deficit decays to |U N (x)|/2 is denoted as r u , which can be chosen as the characteristic length in the spanwise direction.The self-similarity of the mean velocity is manifested after normalization by the characteristic length r u (x) and velocity difference |U N (x)|.The maxima of velocity deficit are located around the centerline stretching from the turbine hub to downstream.As shown in Figure 9a,b, the velocity deficit maximum U N (x) decays with the streamwise distance x while the half width r u grows with x.If the velocity deficit U − U ∞ and spanwise distance y are normalized by the velocity deficit maximum |U N |(x) and half width r u , respectively, the spanwise distributions of the mean velocity deficit collapse, indicating the occurrence of a self-similarity.
The streamwise development of the self-similarity is categorized into three stages based on the characteristic scales and the normalized velocity profiles.The thresholds of these stages can be determined from the changes of the growth rate of the half width r u .The thresholds can also be confirmed by the decay rate of the velocity deficit maximum |U N |.Three stages are separated by the dotted lines in Figure 9a,b.Stage I in x/D < 2 is the near wake region, where the velocity deficit profile does not fit a Gaussian distribution, meaning that the flow has not entered the regime of self-similarity yet.The normalized velocity profiles of stage II at 2 < x/D < 5 and stage III at x/D > 5 fit the same Gaussian distribution in the spanwise direction, as shown in Figure 9c.Therefore, stages II and III are in the regime of self-similarity for the mean velocity.However, the change of growth rate at x/D ∼ 5 indicates that an incident happens there.In the literature review of Xie and Archer [30], another three-stage categorization was mentioned: the near wake, the transition region, and the far wake.The locations of the boundaries of the stages are similar between these two categorizations.To clarify, our categorization is based on the quantitative analysis of the characteristic scales.The second stage of our categorization is not a "transition stage" as the self-similarity is observed.The mechanisms behind the establishment and development of the self-similarity are explained in Sections 4.4 and 4.5.

Self-Similarity in Reynolds Stress
The Reynolds stress u i u j plays a vital role in the evolution of mean velocity field in the turbulent wake flow.With a similar process of defining the characteristic scales and normalizing afterwards, the self-similarities are also found for the Reynolds stress components u u and u v in the horizontal plane at the hub height.
For the Reynolds stress component u u , the general shape of its spanwise distribution in the stages II and III is shown in the normalized profile in Figure 10b.It has two peaks at y/D ∼ ±0.5 and a trough at y/D ∼ 0. The characteristic Reynolds stress scale (u u ) N is defined as the difference of u u between the peak and the trough, while the characteristic width scale r u u is defined as half of the spanwise distances between the two peaks.The normalized profiles of u u in the region 4 ≤ x/D ≤ 8 well collapse to the same curve, as shown in Figure 10c.For the region surrounding but outside 4 ≤ x/D ≤ 8, the normalized profiles are in a similar shape but with some deviations.A subtraction of the spanwise averaging between the two peaks, u u = [ 2πru u dr]/[πr 2 ], is necessary in the normalization to achieve the collapse.The normalized Reynolds stress (u u − u u )/(u u ) N in the side y/D > 0 is a bit larger than that in the side y/D < 0. For the streamwise development of the characteristic scales, as shown in Figure 11, a change of the tendency is observed at the boundary of the stages II and III.The length scale r u u , namely the spanwise peak-to-peak distance, stays at a flat level after entering the stage III.
For the Reynolds stress component u v , its general shape of the spanwise profile is shown in Figure 12b.Two different characteristic width scales are found inside and outside the peak/trough pair, respectively, by which an inner region and an outer region are separated.For the inner region, the characteristic Reynolds stress (u v ) N and the characteristic width r u u ,inner are defined as the spanwise distance and the Reynolds stress difference, respectively, between the pair of peak and trough.For the outer region, the characteristic width r u u ,outer is defined as the distance to the peak/trough where the Reynolds stress decays to a half.As shown in Figure 12c, the normalized profiles of u v in the region 3 ≤ x/D ≤ 10 well collapse to the same curve.Derived from the boundary-layer equation shown in Equation ( 23), the normalized profile of u v has an analytical solution where η u v is the normalized spanwise coordinate defined as below, The streamwise development of the characteristic scales is shown in Figure 13.The entry of the stage III causes a noticeable change of tendency only in the inner width scale.
The spanwise gradient of the mean velocity is assumed to be the source of the turbulent kinetic energy (TKE) in the turbine wake.Therefore, the mean velocity deficit scale |U N | is used in the normalization of the Reynolds stress scales in Figures 11a,b  Figure 10.The normalization for the spanwise profile of the Reynolds stress component u u .In (a), the spanwise profiles are normalized by a fixed scale for all downstream distances.In (b), the definitions of the characteristic scales in u u are annotated.In (c), the spanwise profiles in 4 ≤ x/D ≤ 8 are normalized by the characteristic scales.In the normalization, u u is subtracted, which is the average of u u for grid points between the two positive peaks.−0.005 0.000 0.005

Mean Momentum Equation and Boundary-Layer Equation
With certain assumptions and simplifications to the Navier-Stokes equations, the boundary-layer equations provide theoretical solutions for boundary-layer flows.The self-similarities occur after a development stage, which is a long distance from the static axisymmetric body.Johanasson and George [50] reported that in a disk wake flow of Reynolds number Re D = U c D/ν = 26, 400, the self-similarity starts from x/D = 10 for the mean velocity and x/D = 30 for the turbulent intensity.Uberoi and Freymuth [51] reported that in a sphere wake flow of Reynolds number Re D = U c D/ν = 8600, the self-similarity starts from x/D = 50.In our simulation, the self-similarity of streamwise velocity starts from x/D = 2 and the self-similarity of Reynolds stress starts from x/D = 3 at a Reynolds number of Re D = 13, 100, which has a comparable Reynolds number but at a shorter distance.Although the earlier occurrence of the self-similarity in the turbine wake is a bit surprising at first glance, it can be explained by a budget analysis of the mean momentum equation.
In the cylindrical coordinate system (x, r, θ), the continuity equation is The mean momentum equation on cylindrical coordinates is The terms on the left-hand side of Equation ( 22) are the mean flow convections.The last three terms on the right-hand side are effects of the Reynolds stress.In axisymmetric wakes, firstly the terms containing derivatives with respect to angle θ diminish due to the symmetry.Secondly, in turbulent flows, the terms containing molecular viscosity ν are neglected as they are smaller than the terms containing the turbulent viscosity ν t .Thirdly, the pressure gradient −∂p/∂x is assumed to be small in the self-similarity regime that is far away from the body.Lastly, −∂u x u x /∂x, the streamwise derivative of Reynolds stress, is neglected because it is much smaller than the radial derivative term −[∂(ru x u r )/∂r]/r.
Based on the considerations above, a simplified mean momentum equation is obtained as To close the equations, the turbulent viscosity and mixing-length model are used to estimate the Reynolds stress term in Equation ( 23) [2], where ν t is the turbulent viscosity for the axisymmetric wake and α is a non-dimensional constant for the turbulence.The classical solution to the boundary-layer Equations ( 21) and ( 23) where C D is the drag coefficient of the axisymmetric body, λ is the dimensionless expansion rate for the wake width ∆, and x 0 is the virtual starting location of the self-similarity.The relation between r u and ∆ is r u /∆ = √ ln 2. The classical solution assumes that where C ss is a x-independent parameter.As a result, the boundary-layer Equation ( 23) is simplified to an equation of the dimensionless coordinate η, and the self-similarity is obtained.The solution of the mean velocity predicts that Equation (30) indicates that the boundary-layer Equation ( 23) yields a Gaussian distribution of the mean streamwise velocity in the radial direction.
The derivation of the Gaussian distribution for the streamwise mean velocity helps us understand the self-similarity shown in Figure 9c.However, considering the planar boundary layer and the rotation effect, the turbine wake is not a purely axisymmetric flow.The assumptions in the simplification of the momentum equation should be verified for the turbine wake flow, which can be done by examining the mean momentum budgets of Equation ( 22) using the LES data.A temporal averaging is firstly performed at each grid point.Next an azimuthal averaging is performed to acquire the temporal-radial mean field.Figure 15a shows the radial distribution of each mean momentum budget term at x/D = 3.5, which is the center of stage II of the self-similarity.By further applying a radial averaging with the radius as the weighting function, the temporal-azimuthal-radial averaged momentum budget is acquired.As shown in Figure 15b, the dominant terms in the budget of mean momentum are the radial turbulent production term −[∂(ru x u r )/∂r]/r and streamwise convective destruction term −u x ∂u x /∂x.The radial convective production term −u r ∂u x /∂r is the third dominant term, which is smaller than the other two leading terms in the region 2 < x/D < 6 and becomes negligible in the region x/D > 6.The fourth dominant term is the streamwise turbulent destruction term −∂u x u x /∂x, which contributes a small fraction of mean momentum in the region x/D < 5.The pressure gradient term −∂p/∂x is important in the region x/D < 2 and becomes negligible in the region x/D > 3.By overlooking the development of the mean momentum budget terms in the streamwise direction as shown in Figure 15b, it could be found that the leading terms in stage II and III are the streamwise convection term −u x ∂u x /∂x, the radial derivative of Reynolds stress term −[∂(ru x u r )/∂r]/r, and the radial convection term −u r ∂u x /∂r.
The leading terms in the mean momentum Equation ( 22) can also be confirmed by simple calculations based on the self-similarity and the characteristic scales analyzed in Sections 4.2 and 4.3.For the streamwise convection term −u x ∂u x /∂x in the region 4 ≤ x/D ≤ 6, the mean velocity u x is estimated as 0.66U ∞ , which is the difference of the inflow velocity U ∞ and the averaged velocity deficit scale |U N | = 0.34U ∞ .The derivative ∂u x /∂x is estimated as d|U N |/dx ∼ 0.1U ∞ /(2D), with the gradient read from the slope at x/D = 4 and x/D = 6 in Figure 14a.The estimation process is reorganized as below, Other three leading terms found in Figure 15 can be estimated in a similar approximation, as shown below, Here, these four terms are estimated using data in the region 4 ≤ x/D ≤ 6.In the estimation of the third term −u r ∂u x /∂r, v N is the scale of the radial mean velocity u r .The scale v N (x) is defined as the maximum radial mean velocity at a downstream distance x.In the region 3 ≤ x/D ≤ 10, v N (x) changes from −0.009U ∞ to −0.006U ∞ as x increases.As the radial velocity scale is one order smaller than other velocity-like scales, its development along x contains a lot of noises.Thus, it is not shown as other variables of interests.
Compared to the cylindrical averaged results in Figure 15, the characteristic scale-based estimation above provides the same sequence of importance for the terms in the mean momentum budget, The fourth term −∂u x u x /∂x is one order smaller than the third term.The third term −u r ∂u x /∂r is nearly one order smaller than the leading two terms, which makes the leading two terms seem to balance with each other.By ignoring the small terms, three leading terms constitute to the boundary-layer Equation (23), which has a Gaussian solution (30) based on past studies of classical axisymmetric wake.Therefore, it is not surprising to observe that the self-similarity occurs as early as x/D = 2 behind the turbine, compared to the starting point x/D > 50 in classical axisymmetric wake at a similar Reynolds number.
The mean momentum budget analysis also throws light on how to get more accurate theoretical solutions for the mean flow of turbine wake.At stage I (x/D < 2), the pressure gradient and streamwise convection are the dominant terms.At stage II (2 < x/D < 5), the consideration of streamwise turbulent destruction in the boundary-layer equation should give a better solution.In region x/D > 6 of stage III, the boundary-layer Equation ( 23) could be further simplified to a two-term equation where the radial convection is neglected.

The Collision Hypothesis of Wake Meandering and Tip Vortices
Under the perturbations from the turbulent flow upwind, the wake flow of a turbine has a transverse meandering with strong interactions between the wake flow and the ambient flow.The velocity deficit centerline is often used to represent the motion of wake meandering [12,31].
At each streamwise location, the spanwise coordinate of the wake centerline is identified with the velocity minimum in the transversal cross-section.A spanwise averaging window is applied in the search of the velocity minimum to reduce fluctuations.As shown in Figure 16a, the centerline originates from the turbine hub and extends downstream with a gradually increasing spanwise amplitude.By superposing the wakelines at various times in the simulation, the region where the low-speed wake sweeps is acquired, as shown in Figure 16b.Some extreme samples of the wake meandering are filtered out, which leaves the coverage map with a pair of smooth edges.In our simulation, the expansion angle of the wakeline is 11.9 • .Tip vortices are another source of turbulent kinetic energy in the wake flow.They are convected downstream under the influences of the ambient flow in the outer region and the velocity deficit in the inner region.Howard et al. [12] assumes that tip vortices move along y/D = ±0.5, which corresponds to a width of D between tip vortices in a horizontal slice crossing the hub.The spanwise velocity around the tip vortex is simplified as zero in the assumption of [12].However, as the wake expands downstream, the spanwise velocity should point to the central axis to entrain the ambient flow.Considering the continuity Equation ( 21) and the decreasing gradient of the streamwise velocity along x, the amplitude of spanwise velocity around the tip vortex should also decrease with x.Consequently, the tip vortices may have a small inward displacement due to the spanwise velocity of the flow.The inward displacement may have a limit as the amplitude of spanwise velocity decays with x.Hu et al. [52] presented the evolution of the unsteady vortex in turbine wake with phase-locked PIV measurements.Although they did not state explicitly, the inward displacement of tip vortices could be observed in the scattering of the wandering tip vortices.Based on the analyses above, the width of tip vortices may be smaller than the rotor diameter D. By assuming the tip vortices are located at the positive and negative peaks of the vertical vorticity ω z in the horizontal plane at the hub height, the width of tip vortices could be estimated from the mean vorticity field in the simulation.As shown in Figure 17, the width W in region 2 < x/D < 7 is at a steady value 0.8D.
Given the expansion rate of the velocity profiles and the nearly constant width W of the tip vortices, the collision of these two flow structures is expected to occur at downstream location.As shown in Figure 16b, the collision location is estimated to be x/D ∼ 4.4 for W = 0.8D or x/D ∼ 5.4 for W = 1.0D.
The development of TKE also discloses the interactions in the wake flow of turbine.As shown in Figure 18, strong TKE k is observed at x/D ∼ 5, around the estimated location of the collision shown above.The collision leads to generation of TKE.The Reynolds stress τ t = −u x u r is enhanced due to the collision.This further leads to a larger turbulent viscosity ν t in Equation ( 24), a larger dimensionless turbulence parameter α in Equation ( 25), and a larger λ in Equation (28).Noted that λ is the dimensionless expansion rate of the wake width ∆ in Equation ( 27), there is thus a larger growth rate of half wake width r u in Figure 9.As the collision position is located at the boundary of stages II and III of the self-similarity, x/D = 5, the interaction between the meandering velocity deficit and the tip vortices leads to the transition from stage II to stage III.The asymmetry is observed in the distribution of TKE.In two previous experiments, the asymmetry was also observed in either the turbulent intensity [44] or the wake expansion angle [12].In those experiments, the turbine GWS/EP-6030 rotates along −x direction and the activity of interest is stronger in −y.In our simulation, we set the turbine to rotate along +x direction and the activity of interest is stronger in +y.Therefore, the asymmetry might be caused by the rotation of the turbine blades.In Figures 10 and 12, the normalized profile of u u shows an asymmetry while u v does not.

Conclusions
In this study, the wake flow of a model-scale wind turbine is investigated by LES.A hybrid derivative approach sPSMFDM is proposed to alleviate the Gibbs phenomenon in the simulation using the spectral method and the actuator line model.The self-similarity of mean velocity in the turbine wake and its mechanism are investigated.The following conclusions are obtained: The Gibbs phenomenon can be alleviated by the sPSMFDM approach, which smooths the input of the PSM and rectifies the output by the FDM.The accuracy of the simulation is validated through comparisons of flow statistics to wind-tunnel measurements.
By defining the maximum velocity deficit and the half velocity profile width as the characteristic scales of velocity and length, respectively, the self-similarity of mean velocity is illustrated.Based on the development of the self-similarity, a new three-stage categorization of the turbine wake is proposed.Stage I is the developing near wake.The self-similarity starts in stage II where the normalized velocity profiles collapse to the same Gaussian distribution in the spanwise direction.Stage III also obtains the self-similarity, while the expansion rate of wake becomes larger than at stage II.The self-similarities of the Reynolds stresses u u and u v in the turbine wake are investigated for the first time.The normalized spanwise profiles of the Reynolds stresses collapse to the same curves after defining the characteristic scales properly.The theoretical formula for the normalized u v is provided.
The mechanisms in the establishment and the development of the mean velocity self-similarity are discussed.A mean moment budget analysis is performed using the cylindrical averaging and the characteristic scale-based estimation.The budget analysis shows that, the leading terms affecting the mean momentum are the streamwise convection and the radial gradient of Reynolds stress.It confirms the assumptions in simplifying the boundary-layer equation for the axisymmetric wake.The Gaussian solution of the simplified BLE is the mathematical basis for the self-similarity in the mean velocity starting at stage II.By recalling the condition in the derivation of self-similarity solution.The increment of the turbine wake expansion rate at the boundary between stage II and III is explained.At that boundary, the collision of the tip vortices and the meandering wake deficit enhances the turbulence.The condition in the derivation of the self-similarity indicates that the enhanced Reynolds stress at the collision leads to the increment of the expansion rate of the turbine wake.
The width of the discontinuous region being smoothed is an important factor affecting the numerical error.In this test, we use a fixed ratio strategy, which determines the smoothing region by expanding the peak-to-trough region with a fixed ratio for each spatial resolution.An optimal ratio is found by iterative attempts to acquire the least error at each resolution.The optimal ratios are 1.67, 0.92, 0.51, and 0.24 for grid number 32, 64, 128, and 256, respectively.The trend of the optimal ratio shows that a narrower smoothing region is needed for a spatial discretization of more grid points.
Several schemes for the PSM are tested to show the effect of sPSMFDM and Fourier smoothing.The hybrid approach of pseudo-spectral method and finite-difference method (sPSMFDM) is generally more accurate than the classical spectral method (marked with PSM), as shown in Figure A2b.The schemes using Fourier smoothing-based dealiasing (marked with FS) are generally more accurate than the classical 2/3 rule-based dealiasing (marked with CD).

Figure 1 .
Figure 1.The estimation of aerodynamic lift force F L and drag force F D on an airfoil section.The upstream fluid velocity is decomposed into a parallel component u and a perpendicular component u ⊥ with respect to the rotation axis.The rotation speed is the product of the angular velocity Ω and the distance r to rotation center.U rel is the relative velocity in coordinate system of the airfoil section.

Figure 2 .
Figure2.Illustration of the hybrid approach sPSMFDM.The velocity field at t = 0.159 for an initial value problem of 1D Burgers equation is shown in (a).The original velocity field is denoted by -, the smoothed field is denoted by , and the points for cubic interpolation are denoted by .The spatial derivative ∂u/∂x is plotted in (b).The analytic solution is denoted by -, the result of the hybrid approach sPSMFDM is denoted by +.The result of PSM on a smoothed input (sPSM) is denoted by .The result of FDM for correction of sPSM is denoted by •.The result of PSM on the original velocity field is denoted by × for comparison.

Figure 3 .
Figure 3. Illustration of the computational domain.The wind blows in the x direction.A turbine model is deployed at a distance of 0.75 m from the inlet.A turbulent inflow database is imported in the velocity relaxation zone upstream.
Chord length and twist angle

Figure 5 .
Figure 5.The grid discretization for fluid and turbine in y − z plane.The rectangular cells are for fluids.The red lines inside the turbine geometry are actuator line elements.

Figure 6 .
Figure 6.Comparison of the contours for − ū∂k/∂x, the streamwise advection of turbulent kinetic energy (TKE), calculated by (a) the standard pseudo-spectral method and (b) the proposed sPSMFDM, respectively.

= 5 Figure 7 .
Figure 7. Vertical profiles of time-averaged velocity U(z), turbulence intensity u rms , and time-averaged Reynolds stress −u w at two different distances downstream of the turbine: subplots (a-c) are at x/D = 2 and subplots (d-f) are at x/D = 5.The results from our LES are denoted by -; the measurement data in wind-tunnel experiment of [44] are denoted by •.In (a,d), the inflow profiles are plotted for reference: the LES results of inflow are denoted by −−; the wind-tunnel data of upstream flow are denoted by .The heights of rotor top, rotor center, and rotor bottom are denoted by • • • .

Figure 8 .
Figure 8.The mean velocity of the turbine wake flow.The bold line in (a) represents the rotation region of turbine blades.The velocity along the dashed line in (a) is plotted in (b).

Figure 9 .
Figure 9.The self-similarity of the mean velocity deficit.The characteristic scales of velocity deficit and wake width are shown in (a,b), respectively, where the dotted lines are the boundaries of the three stages of wake development.The normalized spanwise profile of mean velocity deficit at the stages II and III are plotted with circular points in (c) and the dashed line is a Gaussian distribution fitting for the normalized profile.
and 13c.The |U N |-normalized Reynolds stress scales show nearly linear dependence with the streamwise distance, which supports the assumption of TKE source.As |U N | is decreasing with x, the absolute trends of the Reynolds stress scales are more complex than a linear relation.The comparison of four velocity-like scales |U N |, u u , (u u ) N and (u v ) N are shown in Figure 14a,b.By normalization with a x-irrelevant parameter U ∞ , the Reynolds stress scales grow at the starting stage and decay afterwards.The comparison of three length scales r u , r u u , and r u v are shown in Figure 14c.In stage III, the widths of u u and u v are at a stable level, while the width of velocity deficit keeps growing.

Figure 11 .
Figure 11.The streamwise development of the characteristic scales of the Reynolds stress component u u .The spanwise averaging u u is shown in (a).The characteristic Reynolds stress (u v ) N is shown in (b).In (a,b), the Reynolds stress is normalized by the characteristic velocity deficit |U N |.The characteristic width r u u is shown in (c).The dotted lines denote the boundary of the stage II and III of the wake development.

Figure 12 .
Figure 12.The normalization for the spanwise profile of the Reynolds stress component u v .In (a), the spanwise profiles are normalized by a fixed scale for all downstream distances.In (b), the definitions of the characteristic scales in u v are annotated.In (c), the spanwise profiles in 3 ≤ x/D ≤ 10 are normalized by the characteristic scales.The vertical coordinate η u v in (c) is defined by Equation (20).The dashed line in (c) is the theoretical solution Equation(19) for the normalized profile of u v .The dotted lines in (c) are the boundaries of the inner and the outer regions for the spanwise distribution.

Figure 13 .
Figure 13.The streamwise development of the characteristic scales of the Reynolds stress component u v .The characteristic widths of the inner region and outer region are shown in (a,b), respectively, where the dotted line is the boundary of the stage II and III of the wake development.In (c), the characteristic Reynolds stress (u v ) N is normalized by the characteristic velocity deficit |U N |.

Figure 14 .
Figure 14.The of the characteristic scales.Four velocity-like scales are shown in (a), three of which are shown in (b) as an enlarged view.The length scales are shown in (c).The scales for the mean velocity deficit U − U ∞ , Reynolds stress components u u , and u v are denoted by -, +, and •, respectively.The square root of the spanwise mean of the Reynolds stress u u is denoted by -.

Figure 16 .
Figure 16.Wake meandering.The instantaneous centerline of the meandering velocity deficit is shown in (a).The superposition of centerlines is shown in (b).The dotted and dashed lines in (b) are the edges extracted from the superposition of centerlines at all timesteps, while the circles are the wakeline at 25 random timesteps.The dashed lines A and B in (b) are the locations where the centerline may collide with the tip vortices.

Figure 17 .
Figure 17.The averaged vertical vorticity in the horizontal plane at the hub height.The positive and negative peaks of the vorticity are marked with the dashed lines.

Figure 18 .
Figure 18.The averaged TKE in the horizontal plane at the hub height.

Figure A2 .
Figure A2.The dependence of the averaged error on the grid number in (a) FDM schemes and (b) spectral schemes.The finite-difference schemes are compared in (a), where FDM means the finite-difference method, O followed by a number means the order of accuracy, UP means upwind scheme, CT means central difference, and AS means adaptive stencil scheme.The spectral schemes are compared in (b), where PSM means pseudo-spectral method, sPSMFDM means the hybrid approach of smoothed PSM and FDM, CD means classical 2/3 dealiasing, FS means Fourier smoothing, and ND means no dealiasing is applied.