LES Study of Wake Meandering in Different Atmospheric Stabilities and Its Effects on Wind Turbine Aerodynamics

Abstract

Wake meandering disturbs the stability of the far wake field and thus increases the fatigue loads of downstream wind turbines. A deep understanding of this phenomenon under atmospheric boundary layers and its relation to the structural loads helps to better model the dynamic wake and alleviate adverse effects. A large eddy simulation and an actuator line model are introduced in the present work to simulate the wake field and aerodynamic loads of wind turbines with different longitudinal spacings. By temporal filtering and the gaussian fitting method, the wake center and edge are precisely defined, and the dynamic wake characteristics, including the wake width, oscillation amplitude, and frequency, are described based on the statistical data of the simulated flow field. Results reveal that the wake meandering is caused by both large-scale atmospheric structure and the unstable vortex shed from the rotor because two distinct meandering frequency ranges are detected. As the atmosphere instability increases, the former becomes the dominant inducing factor of the meandering movements. Further, the analysis of the correlation between the inflow characteristics and the wake deflection shows that the Taylor hypothesis remains valid within a distance of over a thousand meters under both neutral and convective boundary layers, proving the feasibility of using this hypothesis for wake evolution prediction. In addition, our study shows that the fluctuation of blade root moment and yaw moment is significantly intensified by the meandering wake, with their standard deviation is augmenting by over two times under both atmospheric conditions. The power spectrum illustrates that the component with rotor rotation frequency of the former is sensible to the wake effect, but for the latter, the power spectrum density of all frequencies is increased under the meandering wake. These indicate that the fatigue loads will be underestimated without considering the wake meandering effect. Moreover, the high correlation between the wake deflection and yaw moment implies that we can predict yaw moment based on the incoming flow information with high accuracy.

1. Introduction

The wind energy industry has seen rapid progress in the last decade, contributing a global cumulative power output of 591GW [1]. The intensification of climate change and the development of related technologies will lead to a continuous prosperity of wind power in the foreseeable future. A horizontal axis wind turbine (HAWT) is the most common wind energy converter in commercial wind farms, and every working HAWT will inevitably leave a tube-shaped wake region characterized by relative low wind speed and high turbulence intensity, which extends over a 10-diameter distance and reduces the performance of turbines located downstream [2,3]. Some analytical wake models [4,5,6] have been proposed to describe the velocity distribution behind wind turbines, aiming to predict and optimize the output of wind farms. Although widely applied to engineering practice over past decade, these models only provide static wake field results at a large time scale and thus cannot support a real-time cooperated wind farm control strategy. Therefore, further research is required to gain a deep understanding of the dynamic characteristics of wind turbine wake.

Wake meandering—large-scale movement of the whole wake flow in both horizontal and vertical directions—usually appears after the full development of a jet or wake flow, such as the large-scale structure in the two-dimension jet flow [7,8], the fin wake of a submarine [9]. and especially the wake of hydro turbines or wind turbines under turbulent ambient flow [10,11,12]. In the wind farm environment, the wiggling turbine wake deteriorates the inflow conditions of downstream rotors by partly shadowing them with varying degrees and hence remarkably increases the output instability and fatigue loads of wind turbines [13]. Numerous wind-tunnel and numerical experiments have been conducted to investigate this phenomenon, but there are still disagreements on the origin and evolution mechanism of wake meandering. Some studies emphasized the contribution of the unconditional instabilities of the shedding vortex to the meanders in the far wake region [14,15,16]. Kang et al. [17] and Foti et al. [18,19] also indicated that the interaction of the spiral flow structure behind the hub with the rotating direction opposite to the turbine rotor and outer shear layer developed from the shed tip vortex could intensify the meandering movements. Moreover, some experimental studies [11,20] showed a regular Strouhal number (St) between 0.2 and 0.3 for the periodic wiggle of the far wake behind a rotor in various working regimes. However, other researchers pointed out that the relatively idealized inflow condition in the above experiments (e.g., uniform or homogenous turbulent inflow) could favor the evolution of the tip and root vortex and amplify their influence on the behaviors of the far wake [21]. España et al. [22,23,24] observed the far wake of an actuator disk and a rotating turbine at 1:400 geometric scale under quasi-atmospheric boundary layer and isotropic turbulent inflow generated in a wind tunnel, and the movement of the whole wake was only detected under ABL condition. They thus negated the influence of the periodic vortex shedding and considered the turbulent structure with the scale comparable to the wake width as the main reason for the meandering phenomenon, which is also supported by the wake field data gained by the light detection and ranging (LIDAR) system in the work of Bingöl, Trujillo et al. [25,26]. Based on the assumption that the wake movement is led by the lateral and vertical velocity variation of the large-scale turbulent vortex in the incoming flow, some dynamic wake models [27,28,29,30,31] were proposed and improved to describe the wake meandering phenomenon. They are considered to be promising tools for a wind farm optimal control strategy, but certain biases caused by the limit of the Taylor frozen hypothesis and the nonlinear aspects of the relation between upstream velocity and wake position, as indicated by Muller et al. [32], require further analysis.

Furthermore, the stability of the atmospheric boundary layer has a significant impact on the evolution of the wake, because the wind shear, turbulent intensity, integral length scale, and other important inflow properties vary with different atmospheric stratifications [33,34,35,36]. Coudou [37] and Foti [38] produced the dynamic wakes under quasi boundary layer flow generated in a wind tunnel and found a meandering frequency with stable St independent of rotor design. They thus believed that the meandering in the far wake is the amplification of the intrinsic wake instabilities by large scale turbulent eddies within the inflow. Results of the recent theoretical studies seem to further confirm this viewpoint: Mao and Sørensen [39] found the optimal inflow perturbations to induce meandering and illustrated that only perturbations with a certain frequency range could be amplified in far wake. Gupta and Wan [40] developed low-order model based on the amplification mechanism to reproduce the main qualitative features of wake meandering. Nevertheless, these studies were limited within either uniform, sinusoidal inflow, or artificial boundary layer flow generated in the wind tunnel, which cannot fully represent the characteristics of the flow environment in the atmosphere. Actually, the meandering energy in frequency domain is not concentrated but distributed in two main ranges that vary with different atmospheric stabilities. This was preliminarily validated by the measurements of Heisel et al. [41] and will be discussed in detail in the present work. Besides, only a few studies have been conducted to illustrate the variation of the statistical properties of the wake meandering under different ABL conditions and to estimate the resulting extra aerodynamic loads of the wind turbines downstream. Abkar and Porté-Agel [42] discussed the effect of atmospheric stratifications on the mean velocity deficit distribution, turbulence statistics, and meandering intensity in the wake region, but the interaction between the wake and atmospheric flow structures and the main cause of wake meandering were not mentioned. Some recent field studies [43,44] provided valuable wake information under various atmospheric conditions, but dynamic features such as meandering spectrum and wake structure evolution were not available due to the limitation of the measurement equipment.

One of the primary tasks at present for the wind energy community is to realize the optimal control of overall wind farm performance [45], and this is challenged by the complexity of both the dynamic wake itself and its interaction with various ambient conditions [46]. Under such a circumstance, the present work aims to analyze and characterize wake dynamics through statistical data gained from simulations, explain the mechanism behind the wake meandering under different atmospheric stabilities, and identify the spectra and magnitudes of main structural loads caused by this phenomenon. We simulate the wind turbine wakes under neutral and unstable atmospheric boundary layer flows with a low surface roughness length corresponding to the sea-surface condition by using a large eddy simulation solver developed on the OpenFOAM open source library, and the actuator line model is introduced to compute the aerodynamic forces generated by the wind turbines. Detailed information of wake meandering (e.g., amplitude, spectrum, probability distribution) in both horizontal and vertical directions and various downstream positions is presented, and the structural loading instability of the downstream turbines is also investigated through correlation curves and power spectrum analysis.

2. Materials and Methods

2.1. Numerical Methods

2.1.1. Governing Equations

Compared with the Reynolds-Averaged Navier-Stokes (RANS) method, LES shows superiority in dealing with unsteady multi-scale vortex structures and the resulting turbulence in the flow field [47,48] and is competent to reproduce the meandering motions of the wind turbine wake [49]. Through the spatial filtering technique, LES divides the flow field into the resolved part and the subgrid-scale (SGS) part. The flow information of the former is obtained by solving the filtered continuity equation and incompressible Navier-Stokes equation as follows:

$$\frac{\partial {\overline{u}}_i}{\partial x_i}=0$$

(1)

$$\frac{\partial {\overline{u}}_i}{\partial t}+{\overline{u}}_j\frac{\partial {\overline{u}}_i}{\partial x_j}=-\underset{I}{\underbrace{\frac{1}{\rho }\frac{\partial }{\partial x_i}{p}_d(x,y)}}-\underset{II}{\underbrace{\frac{\partial \tilde{p}}{\partial x_i}}}-\underset{III}{\underbrace{\frac{\partial {\tau }_{ij}^D}{\partial x_j}}}-\underset{IV}{\underbrace{2{\epsilon }_{i3k}{\Omega }_3{\overline{u}}_k}}+\underset{V}{\underbrace{g(\frac{\overline{\theta }-{\theta }_0}{{\theta }_0}){\delta }_{i3}}}+\underset{VI}{\underbrace{\frac{1}{\rho }{F}_i}}.$$

(2)

The overbar represents the spatial filtering and thus the resolved-scale velocity is expressed by ${\overline{u}}_i=u_i-u_i^{\prime }$ ($u_i^{\prime }$ is the SGS velocity). Term I is an artificially added background pressure gradient, which has a linear distribution at the horizontal plane and drives the atmospheric boundary layer flow to the prescribed state, where $\rho =1.225{\mathrm{kg}/\mathrm{m}}^3$ is the constant density of the air. $\tilde{p}$ in term II consists of two parts—the resolved pressure subtracting the driving pressure normalized by $\rho$ and one third of the stress tensor trace, i.e., $\tilde{p}=(\overline{p}-p_d(x,y)+\rho gz)/\rho +\frac{1}{3}{\tau }_{kk}$. In term III, ${\tau }_{ij}$ represents the fluid stress tensor, caused by both viscous and subgrid effect, and ${\tau }_{ij}^D={\tau }_{ij}-{\delta }_{ij}{\tau }_{kk}/3$, where ${\delta }_{ij}$ is the Kronecker delta. In order to consider the self-rotation effect of the earth, the Coriolis force is calculated by term IV, in which ${\epsilon }_{i3k}$ the alternating tensor and ${\Omega }_3=\omega [0,\cos (\varphi ),\sin (\varphi )]$ is the rotation rate vector, with the planetary rotation rate $\omega =7.27\times {10}^{-5}\mathrm{rad}/\mathrm{s}$. Since the density $\rho$ is set to be constant in the whole flow field, the buoyancy effect caused by the spatial temperature variation is simulated by the Boussinesq approximation in term V, where $\overline{\theta }$ is the local resolved potential temperature and the reference temperature ${\theta }_0$ is set to 300 K. The body force $F_i$ in the last term of the N-S equation represents the lift force and drag force generated by the turbine blades. The deviatoric part of the fluid stress appeared in term III that reflects the viscosity and the contribution of the small-scale flow should be modeled and in this work the linear-gradient diffusion assumption and the Smargorinsky eddy viscosity model [50] are introduced to deal with ${\tau }_{ij}^D$,

$${\tau }_{ij}^D=-2{\upsilon }^{SGS}{\overline{S}}_{ij}$$

(3)

$${\overline{S}}_{ij}=\frac{1}{2}(\frac{\partial {\overline{u}}_i}{\partial x_j}+\frac{\partial {\overline{u}}_j}{\partial x_i}).$$

(4)

The subgrid viscosity ${\upsilon }^{SGS}$ is computed as

$${\upsilon }^{SGS}={(C_s\Delta )}^2{(2{\overline{S}}_{ij}{\overline{S}}_{ij})}^{1/2},$$

(5)

where the filter width $\Delta ={(\Delta x\Delta y\Delta z)}^{1/3}$ and the Smargorinsky coefficient $C_s$ is 0.14. Besides, the transport equation of the potential temperature needs to be solved:

$$\frac{\partial \overline{\theta }}{\partial t}+{\overline{u}}_i\frac{\partial \overline{\theta }}{\partial x_i}=-\frac{\partial q_i}{\partial x_i}.$$

(6)

The temperature source term in the right side is also modeled considering the SGS effect:

$$q_i=-\frac{{\upsilon }^{SGS}}{{\mathrm{Pr}}_t}\frac{\partial \overline{\theta }}{\partial x_i},$$

(7)

where the turbulent Prandtl number ${\mathrm{Pr}}_t$ is set to 1/3 in the present work. It should be mentioned that using subgrid-scale models [51,52] that compute the model coefficient locally and dynamically could be better choice, but some comparative numerical studies [53,54,55] found that when having proper mesh and actuator line resolution, different SGS models show small impact on the wake flow statistics and the performance of standard Smagorinsky SGS model can be effectively improved by choosing the Smargorinsky coefficient appropriately. Therefore, the $C_s$ and ${\mathrm{Pr}}_t$ value close to those in similar work [56,57] is adopted to obtain enough accuracy while economizing on the computational cost. In addition, ${\mathrm{Pr}}_t$ actually has a positive correlation with the local flow stability. For simulation of stably stratified boundary layer where flow dynamics are more sensitive to the subgrid-scale motions, more sophisticated subgrid (SGS) model [58] will be necessary, which is the reason why only the wind fields in neutral boundary layer (NBL) and convective boundary layer (CBL) are simulated in the present study.

Based on the open source C++ libraries OpenFOAM of version 2.3.1, the govern equations are discretized on unstructured and collocated meshes with finite volume method. Unlike all other variables stored at cell centers, the subgrid quantities are directly computed on cell faces. In this way the unphysical friction near the ground can be mitigated [56]. The quantities on the cell faces required in the divergence calculation are obtained through linear interpolation from neighbor cells. The oscillated pressure field on the non-staggered grids is avoided by the velocity flux interpolation technique [59]. The pressure implicit splitting operation (PISO) method [60], with second-order backward time discretization, is introduced to realize time advancement, which consists of velocity prediction and correction based on the old pressure field and continuity rule. Temperature, SGS quantities, and body forces are updated after the end of the PISO algorithm.

2.1.2. Actuator Line Model

Since the main purpose of the present work is to study the wake evolution and especially the wake meandering phenomenon in the far wake, the enormous computation due to the high resolution near the surface of turbine blades is not necessary and can be avoided by using the actuator line model (ALM) proposed by Sørensen and Shen [61] instead of a geometric wind turbine. By dividing the wind turbine blades into tens of airfoil segments, ALM replaces each segment with the body force calculated according to the local inflow data and the airfoil aerodynamic performance table.

$$U_{rel}=\sqrt{U_x^2+{(\Omega r-U_{\theta })}^2}$$

(8)

$$L=\frac{1}{2}{C}_l\rho U_{rel}{}^{2}cdr$$

(9)

$$D=\frac{1}{2}{C}_d\rho U_{rel}{}^{2}cdr$$

(10)

The geometric relationship among the local velocity, angle of attack and the aerodynamic force is depicted in Figure 1. $U_x$, $U_{\theta }$, $\Omega$ are axial velocity, tangential velocity and rotation speed of the rotor respectively. The lift coefficient $C_l$ and drag coefficient $C_d$ are functions of the angle of attack $\alpha$. Once the relative velocity and aerodynamic coefficients are obtained through iterative calculation, the lift and drag force can be then expressed by Equations (9) and (10), respectively, where $c$ is the chord length and $dr$ is the width of the airfoil segment.

In this way, every actuator segment is represented by the point body force $f$, the vector sum of $L$ and $D$. The concentrated body force should be projected smoothly by a regularization kernel form a point to a sphere-shaped region before being exerted onto the flow field in order to avoid the possible numerical errors. Therefore, the aerodynamic force generated by all the blade segments at $(x,y,z,t)$ is calculated as

$$F(x,y,z,t)=-{\displaystyle \sum _{i=1}^N{f}_i(x_i,y_i,z_i,t)\frac{1}{{\epsilon }^3{\pi }^{3/2}}}\exp [-{(\frac{d_i}{\epsilon })}^2],$$

(11)

in which $N$ is the total number of the actuator segments and $i$ is the segment index. $d_i$ is the distance between $(x,y,z)$ and the location of the $i\mathrm{th}$ actuator point. $\epsilon$ is a constant parameter deciding the width of the projection region, which has an important influence on the computational quality. In the present work, $\epsilon$ is set equal to 4.0, twice of the local grid side length, which is the recommended value that insures a stable numerical solution [62].

2.2. Simulation Setup

In order to generate the atmospheric boundary layer inflow, the precursor-successor method is introduced, i.e., simulating the whole evolution progress from an initial laminar shear flow to the fully developed turbulent boundary layer state in a precursor case and then collecting the section flow information as the ABL inflow data of the wind farm simulation (successor case). In this way, the coherent large-scale turbulent structures in the ABL, which have significant influence on the wake meandering, can be effectually reproduced and the statistical features of the generated turbulence keep stable from the inlet plane to the end of the domain, guaranteeing a good atmospheric self-balance property.

2.2.1. Precursor Simulation

When deciding the dimensions of the computational domain, we should ensure that the atmospheric turbulent structures with scales ranging from several meters to over a kilometer can be generated. Besides, the streamwise length of the domain should be long enough so that the far wake behaviors of wind turbines can be simulated and studied. Therefore, the length, width and height of the domain are respectively set to 2016 m, 1008 m, and 1008 m. The whole computational domain is meshed uniformly in all three dimensions with a 8 m resolution, i.e., a grid number of $252\times 126\times 126$ in x, y and z axis, which will also be the background mesh of the successor case.

The aim of the precursor case is to simulate two atmospheric boundary layers with different stabilities (NBL and CBL). The evolution from the initial laminar flow to a fully developed ABL flow actually needs a very large streamwise distance, which is far beyond the present domain length. The solution is to apply the cyclic conditions to the four vertical boundary planes. Besides, the upper boundary is subject to the slip condition, meaning no flux and vertical velocity gradient in this plane. The direct simulation of the flow over a rough ground with high Reynolds number is computationally unacceptable, so in the lower boundary, the Moeng’s wall model [63] is introduced to estimate the surface stress on the basis of the velocity near the ground and the given surface roughness length $z_0=0.001\mathrm{m}$, corresponding to the situation over a calm sea surface. The mean wind speed at the hub height plane is the rated wind speed of the wind turbine $u_{hub}=11.4\mathrm{m}/\mathrm{s}$. It should be mentioned that if the main flow direction is aligned with the x-axis, the turbulent structures under cyclic boundary condition could be, somehow, laterally restrained in a limited region, leading to a non-uniform distribution of the average velocity in the horizontal plane, so the wind direction is set at 15.52° counterclockwise angle to the x-axis to prevent this phenomenon and meanwhile maintain an enough streamwise length in the domain. In addition, the surface temperature flux is prescribed to $0.00\mathrm{K}\cdot \mathrm{m}/\mathrm{s}$ and $-0.04\mathrm{K}\cdot \mathrm{m}/\mathrm{s}$, respectively, in the neutral and convective cases. From 0 m to 900 m in the z-direction, the potential temperature is 300 K and then linearly increases to 308 K in the next 100 m distance upward, forming a capping inversion region at the top of the boundary layer.

2.2.2. Successor Simulation

With the background mesh being the same as the Cartesian grids in precursor case, two refinement processes are imposed to the successor computational domain. According to the mesh resolution, the domain is divided into three regions, marked I, II, and III. Starting from $1D$ ($D$ represents the diameter of the turbine rotor) ahead of the first turbine to its $13D$ downstream position in the wind direction and extending $1.5D$ distance outward based on the hub center in lateral and vertical directions, Region II covers the wake flow of wind turbines with a mesh resolution of 4 m. Region III has the finest 2-m resolution, including the vicinity of each turbine rotor. The details of the wind farm layout and the scales of the computational domain are shown in Figure 2.

All the wind turbine used in the present study are the NREL-5MW baseline wind turbines [64], of which the gross properties are listed below in

Table 1. Gross properties of NREL-5MW.

Item	Value	Unit
Rating	5.0	MW
Rotor orientation	Upwind	-
Blade number	3	-
Rotor diameter	126	m
Hub height	90	m
Rated wind speed	11.4	m/s
Rated rotation speed	12.1	rpm

.

Under each ABL inflow condition, four successor cases were setup, among which one contains only a single wind turbine located nearly 3D streamwise distance away from the inlet plane, and the other three cases include two turbines with the second one located $\Delta x=5D,7D,9D$ downstream of the first turbine respectively. For convenience in postprocessing and displaying the flow information, in the rest of the work a new Cartesian coordinate system with its origin based on the location of the wind turbine and x-axis aligned with the wind direction is used and this coordinate system is also shown in Figure 2. Besides,

Table 2. Case number and wind farm layout.

Case	ABL Stability	Number of Turbines	$\Delta \mathit{x}$
1	Neutral	1	-
2	Neutral	2	5D
3	Neutral	2	7D
4	Neutral	2	9D
5	Unstable	1	-
6	Unstable	2	5D
7	Unstable	2	7D
8	Unstable	2	9D

notes the atmospheric environments and the locations of the wind turbines and gives a case number to each successor case in order to avoid confusion in the comparison and analysis of their flow fields.

3. Results

As mentioned in the last section, in this study, two precursor cases were set up to provide the neutral and convective boundary layer inflows for the successor cases. The simulations of both cases first run for 18,000 s with a time step of 0.5 s in order to ensure that the velocity profile in the boundary layer reaches a statistically stable state. Then, the next stage of computation lasts for 840 s, and the time step was changed to 0.02 s, which accords with that in the successor cases. This insures that the displacement of the blade tip point in one time step will not exceed the grid side length, which is important to avoid numerical divergence [62]. The flow information of this stage is stored as the inflow data of the wind farm simulation. Considering that the propagation and full development of the wind turbine wake need at most $12D/u_c\approx 223s$ (assuming the convective velocity of the wake $u_c$ is the minimum time-average wake center velocity 7 m/s), only the flow field data of the last 600 s out of 840 s is used for the analysis in the rest of this work.

3.1. Verification of the Atmospheric Boundary Layer Flow

The time-average velocity profiles and the turbulence intensity profiles in three dimensions of these two types of ABL flow are plotted in Figure 3, noting that the turbulence intensity (TI) is defined as follows:

$$T{I}_i(z)=\frac{\sqrt{\overline{{(U_i(z)-{\overline{U}}_i(z))}^2}}}{\overline{u}(z)}=\frac{{\sigma }_i(z)}{u_0(z)}.$$

(12)

Here, the overbar means time average. $U_i(z)(i=x,y,z)$ is the wind speed in three directions at the height of $z$, and the lower case $u$ represents the streamwise wind speed, which is aligned with x-axis in the new coordinate system. Since the ABL flow reaches the quasi-equilibrium after 18,000 s, the time-average velocity $\overline{u}(z)$ is considered as the ambient wind speed $u_0(z)$ in the wind farm simulation cases.

In Figure 3a, the blue and red solid lines represent respectively the velocity profiles in NBL and CBL. Theoretically, in the neutral situation the velocity distribution in the Prandtl layer, which should be no more than 100 m high, is subject to the following logarithmic law [35]:

$$\frac{\overline{u}(z)}{u_{*}}=\frac{1}{\kappa }\log (\frac{z}{z_0}),$$

(13)

where $\kappa =0.4$ is Karman constant and $z_0=0.001\mathrm{m}$ is the roughness length. Fitting the NBL velocity profile from 0 to 100 m with Equation (13) results in a rational friction velocity $u_{*}=0.398$. The fitted curve is also plotted as black dashed line in Figure 3a, showing that the wind speed profile in the Prandtl layer accords well with the logarithmic law. The use of the wall model instead of direct numerical simulation accounts for the deviation within a few meters above the ground. The hub height velocity in both cases are consistent with the rated wind speed of the wind turbine and the differences between the velocity at the rotor tip and the bottom are 2.00 m/s and 0.95 m/s, respectively, in neutral and unstable situations, presenting a higher wind shear in NBL at the wind turbine working height.

The profiles of turbulence intensity (TI) in three directions are shown in Figure 3b. The streamwise TI profiles in NBL and CBL show differences less than 0.5% at various heights and in both cases the lateral and vertical TI profiles are also quite close to each other. However, under the buoyancy effect caused by surface temperature flux, the turbulence intensities in y and z directions in the convective boundary layer are obviously larger than that in neutral case. This signifies that the transversal turbulent structures in CBL will cause more momentum exchange between the ambient flow and the wake and thus remarkably influence the wake evolution process. Figure 4a,b shows the iso-surfaces of the air mass with the vertical velocity $w=1.0\mathrm{m}/\mathrm{s}$ in the flow field. Compared with the less vertical fluctuations in NBL, the surface temperature flux in the unstable boundary layer leads to the numerous huge air masses with upward tendency, which explains the disparity of $T{I}_z$ in Figure 3b. In the unstable condition, under the effect of buoyancy force and velocity shear, the flow stratification is disturbed, and a turbulent vortex appears and entrains the surrounding air, gradually forming the large-scale flow structures in the atmosphere. This phenomenon is reflected through the velocity contours in Figure 4c,d. The large coherent turbulent structures appear in both cases but in CBL the scale and fluctuation are obviously larger, which will dominate the dynamic wake characteristics in the wind farm.

3.2. Wake Characteristics of a Single Wind Turbine

In this part, the simulation results of the wake of a single turbine under different atmospheric boundary layers (i.e., case 1 and case 5) are analyzed and demonstrated.

3.2.1. Time-Average Wake Field

With part of the kinetic energy converted by the wind turbine, the wind speed decreases sharply behind the rotor. This low speed region is elongated and advected downstream by the ambient flow and expands in transverse direction under the mixing effect of the turbulence.

Figure 5 shows the time-average wake velocity deficit field in the hub height x-y plane and middle vertical x-z plane of both case 1 and case 5. Because the nacelle and tower of turbine are not simulated in the present study, the incoming wind flows past the rotor hub and forms a high-speed tube in the center of the near wake region. In Figure 5a,b, a similar wake evolution process is shown under both neutral and unstable conditions. The shear layer formed behind the blade tips and roots widens as it flows downstream, narrowing the high-speed tube until its disappearance at the onset of far wake region. This dividing point, however, locates at about $3D$ downstream of the rotor in case 5, about one diameter ahead of that in case 1, meaning that the momentum exchange between the wake and ambient flow is more intensive in CBL. Moreover, the velocity deficit profiles in far wake region show perfect self-similarity and axial symmetry, which makes it reasonable to reproduce the profiles by a fitting curve and this technique will be used in the next part to define the wake center and edge. In Figure 5c,d, despite the wind shear to varying degrees in the vertical direction, the wake velocity distributions in different positions downstream are quite similar to those in the hub height plane in both cases. The only difference is that the expansion of the wake is blocked by the ground in far wake region, thus slowing the velocity recovery process near the ground.

To further reveal the influence of different atmospheres on the wake evolution, the dimensionless Reynolds stress contours are plotted in Figure 6. The magnitude of the Reynolds stress firstly experiences an ascent stage when the upper and lower shear layers grow until they meet in the center line, and then tends to a smooth distribution along lateral and vertical directions, signifying the formation of a fully developed turbulent wake flow. It is clear shown from the colors in these contours that the transit of the flow momentum is mostly from the ambient filed to the wake field. The growth of the shear layer is mainly due to the $\overline{u^{\prime }{v}^{\prime }}/{\overline{u}}_{hub}^2$ and $\overline{u^{\prime }{w}^{\prime }}/{\overline{u}}_{hub}^2$ in horizontal plane and vertical plane respectively. From Figure 6, the absolute values of $\overline{u^{\prime }{v}^{\prime }}/{\overline{u}}_{hub}^2$ and $\overline{u^{\prime }{w}^{\prime }}/{\overline{u}}_{hub}^2$ in convective condition is obviously higher than those in neutral condition. A larger Reynolds stress accelerates the development of the wake turbulence and thus the velocity recovery process, which is consistent with the faster expansion and velocity redistribution under unstable boundary layer in Figure 5.

3.2.2. Wake Meandering

Wake meandering, the large-scale lateral oscillation of the wake, is universally observed in the wake field, especially in the wind farm under atmospheric environment. The meandering characteristics could be affected by numerous factors such as the working regime of wind turbine, inflow condition, yaw angle, etc., while the present study focuses on the influence of atmospheric stabilities on the temporal and spatial features of the wake meandering phenomenon.

A filtering technique is first used to screen out the high frequent turbulence while preserving the coherent flow structure information by processing the velocity time-history sequence with a time window filter [18],

$$\widehat{u}=\frac{1}{\tau }{\displaystyle {\int }_{t-\tau /2}^{t+\tau /2}u(t)dt},$$

(14)

where the caret means temporally filtering and $\tau$ is the width of the filter window. It is hard to define the exact boundary between the meandering-scale flow and “totally random” turbulence. Foti [18] proposed $0.5D$ to be the characteristic length for wake meandering, and thus the width $\tau$ should be no larger than $0.5D/u_{hub}=5.5s$ to accurately capture the wiggle of wake center. Additionally, Ott et al. [65] indicated that the “scrambling” eddies with a frequency range of $0.02~0.3\mathrm{Hz}$ have a comparable size with the wake width and contribute a lot to the variation of wake cross section. Therefore, $\tau$ is set to $1/0.3\mathrm{s}\approx 0.6T$ to count in these influences on the wake dynamics. Then, the following three methods for quantitively describing the movements of the wake are discussed:

• Parameter $\mu$ of the best fit gaussian curve to the velocity deficit profiles in the horizontal and vertical planes through the rotor hub;
• The gravity center of the velocity deficit field;
• The point of maximum velocity deficit.

To compare the rationality and accuracy, all three methods are used to estimate the wake center coordinates based on the simulation data. Figure 7 illustrates the results of typical velocity deficit fields in the cross section $8D$ downstream. Wake centers calculated by these methods are marked by an “x,” star, and pentagon, respectively. It should be noted that the calculation should be limited to a proper area to reduce the errors caused by turbulence in the ambient flow. Considering the possible magnitudes of the meandering amplitude and wake area, only the flow information inside the region circled by the dashed line (a circle centered at rotor center with a radii equal to $1.2D$) are adopted in the wake center calculation. Generally, results of three methods differ little in neutral boundary layer as shown in Figure 7a. Nevertheless, when the velocity distribution in the wake is strongly disturbed by the atmospheric turbulence in unstable condition, the result of maximum point method becomes unreliable, because the maximum velocity deficit could probably appear far away from the real wake center in such condition. Figure 7b shows that the maximum point locates at the edge of the wake region while the other two methods seem to give close and rational estimates. Therefore, the gaussian fitting method is adopted for the wake center estimation in the rest of this work.

The original flow data and fitting curves in lateral and vertical directions of these two instantaneous fields are plotted in Figure 8. Despite certain deviations caused by local turbulence, the fitting curves provide good description of the velocity distribution and $\mu \pm 2\sqrt{\ln 2}\sigma$ are defined as the wake edges. With the wake data in every cross section processed in this way, the wake deflections $\delta$ (the distance between the instantaneous wake center and rotor center) in lateral and vertical directions at each moment are obtained. A snapshot of a filtered instantaneous wake flow field is displayed in Figure 9, with the wake center and edge trajectory depicted by the solid and dashed lines respectively. The wakes in both NBL and CBL remain the same size as the rotor and stable until $3D$ position downstream. After that, the wake starts to wiggle around the hub center line in both directions, and the wake width varies randomly as advected downstream.

To investigate the meandering intensity at different downstream positions, 10 min flow data of the simulations obtained with a sampling frequency equal to 50 Hz are filtered and fitted, and the resulting wake movements are analyzed in space and frequency domains. The root mean square (RMS) value of wake deflection explicitly demonstrates the wake meandering amplitude and is computed as follows:

$${\delta }_{rms}(x)=\sqrt{\overline{\delta {(x,t)}^2}},{\delta }_h={\mu }_h,{\delta }_v={\mu }_v-h_{hub},$$

(15)

where ${\delta }_h$ and ${\delta }_v$ signify the wake center deflections in horizontal and vertical directions. Figure 10 shows the wake deflection RMS from $2D$ to $12D$ downstream positions. The results demonstrate that, despite the same turbulence intensity profiles in the y and z dimensions, the lateral wake oscillation is 10–50% larger than those in the vertical direction under both NBL and CBL conditions because the existence of the ground blocks, to some degree, the vertical movements of the wake and atmospheric turbulence. Moreover, it should be noted that the atmospheric instability remarkably intensifies the wake meandering, and especially the oscillation amplitude of vertical wake meandering under convective condition, is even higher than the lateral meandering in the neutral condition.

As advected downstream, the wake of wind turbine uninterruptedly interacts with the ambient turbulent flow, and the latter not only deflects the wake but also changes its shape and area, which has an adverse effect on the performance of downstream wind turbines. In this work, the wake area $A_w$ in each downstream cross section is computed at every moment to have a view of its variations. Because it is difficult to delimit the outer edge of the irregular wake shape, the Monte Carlo method is introduced to estimate the area of the wake. Specifically, in each cross section, we calculate the number of points where the velocity deficit exceeds 0.05 and divide it by the total number of mesh nodes. As the mesh resolution increases, the quotient approaches its asymptotic value, which is considered as the instantaneous wake area. Based on these results, the mean values and probability density distribution (PDD) curves of wake area (nondimensionalized by the rotor area $A_r$) at various positions are illustrated in Figure 11.

As shown in Figure 11, in both NBL and CBL cases the average $A_w/A_r$ increases linearly with the growth of shear layer in near wake region and later slows down its growth rate in far wake region. It is also shown that different atmospheric conditions have little influence on the mean value of the wake area. However, the probability density distributions for two cases show big differences. The PDD curves show an important phenomenon rarely mentioned in previous researches, which is that the wind turbine wake does not simply expand as it flows downstream but the size of its cross-section area varies randomly all the time. Furthermore, the probability distribution of the wake area $A_w$ seems to be subject to a gaussian pattern at the near wake region, but in the far wake, this curve shows a more dispersed distribution instead of an obvious dominant peak. The randomicity of the wake area variation increases as the atmospheric boundary layer flow become more unstable. This dynamic feature is relatively difficult to predict because it is a consequence of the interaction between the turbulence originated from fully developed shear layer and the atmosphere.

Frequency analysis can reveal the driving factors of the large-scale transverse movements of the wake and contribute to the estimation of the aerodynamic performance of downstream wind turbines. Through the Fast Fourier Transform of the 10 min wake deflection time-history sequence, the frequency spectra of wake meandering at $x=5D$, $7D$ and $9D$ in both atmospheric conditions are shown in Figure 12, with the top five peaks marked to illustrate the dominant frequencies. The frequency $f$ in x-axis and amplitude $S_f$ in y-axis of the spectrum are nondimensionalized by $St=fD/u_{hub}$ and $S_f/(DT)$ respectively. A distinct dominant peak corresponding to $St=0.29$ is detected at all three positions in Figure 12a. This dominant frequency of wake meandering is also found in the work of Chamorro [11] and Okulov [20], who attributed it to the instability of the helical vortex structures induced by the rotor blades. This peak is also clear in Figure 12b, but it is different from the results in case 1; under convective condition, more meandering energy is distributed in the region where the Strouhal number is around or lower than 0.1. Since $St=0.1$ is equivalent to a meandering period of 110 s, this low frequent wiggle of the wake could only be driven by the large-scale atmospheric turbulence, which indicates that in neural boundary layer with a low roughness, and the influence of tip vortices on the horizontal wake movements is significant, while under an unstable atmosphere, the atmospheric turbulent structures in the inflow become the dominant driven factors.

In addition, in case 1, the frequency spectra of the meandering in horizontal and vertical direction show obvious difference. The tip vortex–driven meandering still can be detected in Figure 12c, but the amplitudes of $St<0.1$ overwhelm those in the relatively high frequency region, signifying the vortex instability is not the main trigger of vertical meandering even under neutral atmosphere. In contrast, the results of Figure 12b,d show considerable similarity because the wake motions are mainly driven by the same atmospheric inflow. Hence, it is plausible to believe that, as the instability of the ABL increases, the atmospheric turbulence, in place of the rotor-generated vortices, becomes the main driven factor of the wake meandering in both horizontal and vertical planes. Moreover, as seen through comparison between the spectra of different positions, the meandering intensity augments as the wake advected downstream but the frequency of oscillation remains nearly the same. The above results clearly show that the wake meandering is not limited to one dominant frequency. Instead, the oscillation energy is distributed to two main frequency ranges as mentioned in the introduction part, which implying the wake meandering consists of two modes, with one related to the amplified unstable shed vortex, the other induced by the large-scale turbulence in the atmosphere. This explains the divergent opinions in the previous studies about the origin of the meandering phenomenon.

In the article of Larsen et al. [27], the wake motion is presumed to act like a passive tracer of the turbulent inflow structures, which, combined with Taylor frozen hypothesis, leads to a relation between the inflow and the downstream wake position, as follows:

$${\delta }_h(x)=v_c{t}_a;{\delta }_v(x)=w_c{t}_a,t_a=\frac{L}{u_a},$$

(16)

where $v_c$ and $w_c$ are the characteristic inflow velocity in y-axis and z-axis, defined as the average velocity of the rotor-corresponding region in the cross section $3D$ ahead of the wind turbine (because the inflow in this position is not contaminated by the rotor-induced velocity). $L$ is the distance between the inflow cross section and downstream wake cross section. $u_a$ is streamwise advective velocity of the wake and thus $t_a$ is the time consumption of the wake advection between these two cross sections. This theory indicates that the transverse inflow velocity and the wake deflection should be subject to a linear relationship to some degree. However, this is only valid when $L$ is within a distance where the Taylor hypothesis holds. To examine the relationship between the inflow and wake meandering and estimate this effective distance, the correlation function of the time-history characteristic inflow velocity sequence and wake center deflection sequence at various downstream positions is established as

$$Cor{r}_h(\Delta t)=\overline{{\tilde{v}}_c(t){\tilde{\delta }}_h(t+\Delta t)};Cor{r}_v(\Delta t)=\overline{{\tilde{w}}_c(t){\tilde{\delta }}_v(t+\Delta t)}.$$

(17)

The tilde means normalization. In this way, the integration of the normalized curve always equals to one, and the strength of the correlation can be reflected by the comparison between the peak values of the correlation curves in different positions.

The results are displayed in Figure 13. From Figure 13a, under neutral atmosphere, the curves of correlation between characteristic inflow velocity and wake deflection show a distinct peak at $\Delta t=100\mathrm{s}$. This proves that there is strong relation between the transverse inflow velocity and wake deflection in the corresponding direction and according to Equation (16), the time offset $\Delta t=100\mathrm{s}$ could be considered as the advection time of the inflow turbulent structure transported from $x=-3D$ (x-axis position of the inflow cross section) to $x=5D$. As the x-axis position of the wake cross section increases, the time offset $\Delta t$ corresponding to the peak increases because of the extension of the advection distance and the peak value of the correlation remains the same level until $x=9D$. Middle-scale oscillations (with periods of tens of seconds) can be observed in the correlation curves in horizontal direction, reflecting the existence of the vortex induced wake meandering, while in vertical direction, oscillations of similar scale are much weaker, which agrees well with the spectrum results in Figure 12a,c. In the last row of Figure 13a, though the peak in horizontal direction is still distinct, the peak value decreases apparently compared with those in the upstream positions, signifying the weakening of correlation. In Figure 13b, the correlation curves show distinctions from those in neutral condition. Firstly, stronger correlation between the inflow turbulent structures and the wake meandering can be confirmed in both directions according to the higher peak values. Secondly, only large-scale heaves are found in the correlation curves, meaning that the vortex-induced wake motion contributes little to the wake meandering under unstable atmosphere. In addition, the curves of horizontal and vertical directions show similar oscillation patterns, indicating the large-scale inflow structures dominate the wake movements in both directions. Moreover, even at $x=12D$ position, the correlation between the inflow and wake meandering is still apparent. This implies that the Taylor hypothesis holds at a long distance, and thus the wake deflection at a certain turbine position can be predicted with accuracy using the upstream wake flow data.

3.3. Influence of Wake Meandering on the Aerodynamic Loads

Wind turbines in a wind farm could be shadowed to various degrees by the wake upstream due to the wake meandering. Thus, the fatigue loads will be inevitably augmented with the high instability of the inflow. In this part, the results of case 2–4 (NBL) and case 6–8 (CBL) are analyzed and demonstrated, where two wind turbines are placed in tandem with three different longitudinal spacings, $5D$, $7D$, and $9D$. Three types of aerodynamic loads are examined.

• Blade-root out-of-plane bending moment $M_{oop}$: the bending moment in tangential direction at blade root;
• Yaw moment $M_{yaw}$: the moment in vertical direction of the yaw bearing at top of the wind turbine tower;
• Low speed shaft torque $M_{lss}$: the torque experienced by the low speed shaft of wind turbine, proportional to the mechanical power.

To investigate the correlation between wake meandering and the three moments of the downstream wind turbine T2, the correlation function is built as

$$Cor{r}_i(\Delta t)=\overline{{\tilde{\delta }}_h(t){\tilde{M}}_i(t+\Delta t)},i=oop,yaw,lss.$$

(18)

The subscript $i$ represents the three different aerodynamic moments. The wake deflection ${\tilde{\delta }}_h(t)$ at position $2D$ ahead of T2 is used. It should be mentioned that, for simplicity, only the horizontal wake deflection is considered in the correlation function. The resulting correlation curves are shown in Figure 14, from which a distinct peak for $M_{yaw}$ at $\Delta t\approx 23\mathrm{s}$ is clearly seen. The time offset corresponds well to the advection time of the wake along the $2D$ longitudinal distance. This proves that the yaw moment of downstream wind turbine is correlated with the horizontal wake meandering in both atmospheric conditions. When the wake wiggles horizontally, the rotor of T2 is partially shadowed by the wake and the asymmetry of its inflow is accentuated. In this case, the non-uniform x-axial aerodynamic thrust results in very large yaw moment. By contrast, the blade root moment and the low speed shaft moment are more sensitive to the instantaneous aerodynamic environment of every blade, which varies with the azimuthal positions when under inflow with strong wind shear and turbulence. Therefore, the relationship between the wake meandering and these two types of moments is masked by the high frequency fluctuation due to the rotation. Moreover, the correlation value for $M_{yaw}$ under a convective condition is higher than that in neutral condition, indicating the influence of wake meandering on the yaw moment of downstream wind turbine increases as the atmospheric boundary layer becomes unstable.

Standard deviation (STD) reflects the dispersion degree of a set of data. A large standard deviation implies that most of the values have large difference from the average. Therefore, this statistical characteristic value is used to demonstrate the instability of the aerodynamic loads. The STD values of $M_{oop}$, $M_{yaw}$, and $M_{lss}$ for both the upstream turbine (T1) and downstream turbine (T2) are illustrated by the bar charts below in Figure 15. The labels of the horizontal axes $5D$, $7D$, and $9D$ represent the longitudinal spacings of the two turbines in different cases. Figure 15a shows the situation in neutral boundary layer. For the upstream turbine, the STD of $M_{oop}$ is the highest among these three types of moments, because the sensitivity of a single blade to the variation of the inflow is obviously higher than that of a rotor. However, when a wind turbine is placed downstream, the meandering wake of upstream turbines significantly aggravates the asymmetry of its inflow. As a result, the instability of $M_{oop}$ and $M_{yaw}$ experiences a sharp increase. Especially when the spacing equals to $5D$, the STD of $M_{yaw}$ for T2 augments by 2.5 times. Even if the longitudinal spacing changes from $5D$ to $9D$, the standard deviations of $M_{oop}$ and $M_{yaw}$ still remain the same magnitude or only see a slight decline. Different from $M_{oop}$ and $M_{yaw}$, the fluctuation intensity of $M_{lss}$ seems to have little relation to the wake effect. The STD of $M_{lss}$ remains around $300\mathrm{kN}\cdot \mathrm{m}$ no matter where the position is. Unlike $M_{oop}$ and $M_{yaw}$, $M_{lss}$ results from tangential forces experienced by the blades instead of axial forces, with the former much lower than the latter. In addition, the torque of the rotor is an accumulated load from all the blades, so its sensitivity to the wake effect is relatively small compared with the other two types of loads. A similar change tendency of these three moments is also found in the CBL condition. The difference is that the low speed shaft torque of both T1 and T2 increases by more than $100\mathrm{kN}\cdot \mathrm{m}$ and the STD of $M_{oop}$ and $M_{yaw}$ for T2 with a spacing of $5D$ are much larger than those in NBL, but the difference narrows as the spacing enlarges. The more turbulent flow environment and thus more strong wake meandering under convective condition could account for the high STD values of T2-5D in Figure 15b, and the rapid decrease with the growth of spacing is due to the fast recovery of the velocity deficit in the wind turbine wake. The above analysis shows the inhomogeneity of the inflow caused by transverse wake motions remarkably augments the fluctuation of blade root moment and yaw moment. This effect becomes more serious when the atmospheric stability decreases. By enlarging the longitudinal spacing, the fluctuation of yaw moment under unstable condition may be alleviated, but this strategy is unhelpful for reducing the fatigue loads when the atmosphere is relatively stable. Because under such circumstance, the wake velocity deficit recoveries slowly and meanwhile the wake meandering is intensified as convected downstream.

The spectrum analysis could better reveal the relationship between the wake meandering phenomenon and different types of structural loads. The power spectra of the three examined structural moments are plotted in Figure 16 and also the top five peaks of every curve are marked for clarity. In Figure 16a, the peaks of T1 (black pentagons) are mainly distributed in low frequency region, where the Strouhal number is lower than 0.1. However, for the downstream wind turbine, the dominant peaks are concentrated on $St=2.23$, responding to the rotor revolution frequency $f_r=0.20\mathrm{Hz}$. Other peaks of $2f_r$ and $3f_r$ (blade-passage frequency) are also distinct in the spectra of $M_{oop}$. This means that the wake effect under neutral condition mainly increases the fluctuation energy of $M_{oop}$ corresponding to the rotor revolution and blade-passage frequencies. By contrast, in Figure 16b, the power spectrum density (PSD) of T2 at all frequencies is generally higher than that of T1. Figure 16c,d show quite similar energy distributions of $M_{yaw}$ in two atmospheric conditions. In the high-frequency region, only the peak of $3f_r$ is detectable because yaw moment is induced by the ensemble of three blades. Different from the spectra of $M_{oop}$, the dominant peaks of yaw moments agrees well with the meandering frequencies, due to the strong correlation between $M_{yaw}$ and horizontal wiggles of the wake. Again, the shifting of the peak frequencies from $0.1<St<0.3$ in NBL to $St<0.1$ in CBL signifies the weakening of shedding-vortex-driven wake meandering with the decrease of atmospheric stability. Besides, the PSD of T2 is much larger than that for T1 and slowly declines as T2 is placed farther away from T1, which agrees well with the results in Figure 15. For the same reason, in high frequency region we can only see the peak corresponding to $3f_r$ in the spectra of $M_{lss}$. The distribution curves of T1 and T2 in Figure 16e,f are nearly the same in the high frequency region where $St>1$. Besides, in both cases, especially under convective condition, the dominant peaks of T1 are all distributed at the range where $0.01<St<0.1$, indicating the $M_{lss}$ fluctuation of T1 is mainly affected by relative large-scale flow structures. However, for downstream turbines, the peaks move rightwards, of which some frequencies are higher than the meandering frequencies. In fact, the torque of a turbine is related to the total energy flux in the rotor region. Therefore, it is not only sensible to the wake deflection but also the variation of the wake area shown in Figure 11, which accounts for the increase of fluctuation energy of $0.1<St<0.3$.

4. Discussion

The main purpose of the present work is to find out the cause of the wake meandering phenomenon, analyze the variation of its statistical characteristics under different atmospheric stabilities, and demonstrate its effect on the structural loads of a wind turbine. Using the large eddy simulation (LES) method, wind turbine wakes under neutral and convective boundary layer conditions are simulated by the precursor-successor methodology.

From the results presented in part 3.1, it is shown that the wind regimes in neutral and unstable boundary layers have large difference even under the same low surface roughness length $z_0=0.001\mathrm{m}$. Under the effect of surface temperature flux, the convective boundary layer forms a wind regime with less shear and higher turbulence intensity compared with that in neutral condition, leading to the larger Reynolds stress and thus faster velocity recovery in the wind turbine wake flow. The real-time preview control strategy needs the accurate capture of the dynamic wake characteristics to reconstruct to a certain extent the inhomogeneous inflow field of every wind turbine [66,67]. The wake flow structure corresponding to the meandering or larger scale is the most important information and feasible to model [65]. Therefore, the rest of this work is focused on the temporally filtered wake field where only the coherent flow structures are conserved.

It is firstly observed that the unstable atmosphere significantly intensifies the large-scale wake motion in both lateral and vertical directions. It should be noted that the amplitude of vertical meandering in CBL is even larger than the lateral meandering amplitude in NBL. Although this is not consistent with Abkar and Porté-Agel’s study [42], which could be attributed to the differences of the higher freestream wind speed and lower surface roughness in the present work, the result implies the necessity to consider the vertical wake movements in the development of a dynamic wake model. Then, a dynamic wake feature hardly noticed by previous researches is discovered. The cross-section area of wake in a certain downstream distance randomly varies at any time especially in far wake region under unstable atmosphere. This will become the main source of the uncertainty in the power and load prediction. For example, the rotor-effective wind speed, an important quantity for pitch controller [68] defined as the average longitudinal wind speed component over the entire rotor plane, will be obviously disturbed by this phenomenon. In addition, spectrum analysis shows that the meandering phenomenon is comprised of two parts respectively attributed to the shed vortex and the large-scale inflow structures. The former could be dominant under relatively stable ambient flow, but the latter factor becomes more important under convective atmosphere, which provides the reason for the different arguments in related studies [19,22,23,24,25,26,37,38]. Moreover, strong correlation between the inflow characteristic velocity and wake deflection under two kinds of atmospheres is detected. These results, to our knowledge, for the first time investigate the effective distance of the Taylor hypothesis under different atmospheric stabilities and confirm the strong correlation between the turbulent inflow structures and wake meandering. This provides a feasible means to predict the evolution of the wake movements for improving the performance of wind turbine under inhomogeneous inflow—the transverse displacement of a wake could be calculated by the multiplication of a constant characteristic velocity and the advection time.

The last part discusses the relation between the horizontal wake meandering and three kinds of structural moments of the downstream wind turbine. The yaw moment is found to be highly correlated with wake meandering, meaning that it is possible to predict yaw moment according to the upstream inflow data within an appropriate distance. Furthermore, the $M_{yaw}$ sensor would be a perfect indicator of the wake position and thus help to adjust the result of dynamic wake model in control system. Besides, the analysis of STD (see Figure 15) indicates the wake meandering remarkably exacerbates the instability of the blade root moment and the nacelle yaw moment, while the engineering load estimation tools are usually not able to directly reproduce the wake meandering effect. Therefore, the inflow data under various atmospheric stratifications of every downstream turbine gained by nacelle-mounted measurement equipment will be very important information to accurately predict the wake meandering-induced fatigue loads. Finally, the power spectrum density curves further reveal the effect of dynamic wake on different structural moments. For the downstream turbine suffered from wake effect, only the $M_{oop}$ oscillation corresponding to the rotor revolution frequency $f_r$ is augmented, while the power spectral density on the whole frequency domain for $M_{yaw}$ increases obviously. Meanwhile, the dominant peaks of $M_{yaw}$ and $M_{lss}$ are more related to the wake meandering frequencies.

5. Conclusions

In this work, the wake meandering phenomenon under different atmospheric conditions and its effect on the structural loads of a wind turbine were characterized and analyzed through large eddy simulations. We explained the mechanism of wake meandering by spectrum and correlation analysis and meanwhile investigated the frequency and amplitude of wake movements in lateral and vertical directions. The shed vortex-induced meandering is only obvious in the horizontal direction under NBL condition and the large-scale inflow structure becomes the dominant factor in unstable atmospheric condition. Furthermore, this study also confirmed the validity of Taylor hypothesis within at least a distance of one thousand meters in the atmospheric environment, because the peak value of the correlation function between inflow characteristic velocity and wake deflection remains high until $10D$ downstream. Thus, this hypothesis could become a useful tool for wake evolution modeling and wake-affected inflow modeling. In addition, wake meandering–induced load fluctuations were studied, illustrating the remarkable increase of the fatigue loads at the blade root and yaw system caused by inflow inhomogeneity.

These results will contribute to the development of an advanced control-oriented dynamic wake model and the analysis of the structural fatigue characteristics of wind turbines affected by wake effect. Nevertheless, the atmosphere in fact is more changeable and complex than those in the simulations and the layout and working regimes of wind turbines in a real wind farm are not limited to a tandem style and a rated condition. Therefore, more numerical and experimental studies will be conducted to widely learn the wake meandering phenomenon under other atmospheric conditions and find mathematical models and algorithms to predict, based on the Lidar-probed inflow data, the power output, and the fatigue loads of downstream turbines, which will play a significant role in the control and optimization process of intelligent wind farms in the future.