Investigation into Yaw Motion Influence of Horizontal-Axis Wind Turbine on Wake Flow Using LBM-LES ^†

Abstract

The dynamic yaw motion of the wind turbine will affect the overall aerodynamic performance of the impeller and the corresponding wake flow, but the current research on this issue is inadequate. Thus, it is very necessary to study the complicated near-wake aerodynamic behaviors during the yaw process and the closely related blade aerodynamic characteristics. This work utilized the multi-relaxation time lattice Boltzmann (MRT-LBM) model to investigate the integral aerodynamic performance characteristics of the specified impeller and the dynamic changes in the near wake under a sine yawing process, in which the normalized result is adopted to facilitate data comparison and understanding. Moreover, considering the complexity of the wake flows, the large eddy simulation (LES) and wall-adapting local eddy-viscosity (WALE) model are also used in this investigation. The related results indicate that the degree of stability of tip spiral wake in the dynamic yaw condition is inversely related to the absolute value of the change rate of yaw angular speed. When the wind turbine returns to the position with the yaw angle of 0 (deg) around, the linearized migration of tip vortex is changed, and the speed loss in the wake center is reduced at about the normalized velocity of 0.27, and another transverse expansion appeared. The directional inducing downstream of the impeller sweep surface for tip vortex is clearly reflected on the entering side and the exiting side. Additionally, the features of the static pressure on the blade surface and the overall aerodynamic effects of the impeller are also discussed, respectively.

1. Introduction

In the atmospheric environment, wind has a strong random uncertainty, especially in complex terrain areas and urban high-rise residential blocks, there may be complex flow scenarios such as shear wind, side wind, gusts, etc. This will cause the horizontal-axis wind turbine to produce uncertain yaw motion, leading to a series of problems such as reduced energy conversion efficiency, shortened fatigue life, and changes in acoustic radiation characteristics [1,2]. Moreover, the wind speed and wind direction in the local area where the wind turbine is located sometimes change very rapidly, and the regularity of its change is not obvious. At the same time, the yaw adjustment system has inertia lag and other reasons, which causes the horizontal-axis wind turbine to yaw more frequently [3].

Therefore, it is necessary to study the characteristics and changes of wind turbine wake flow field in yaw motion. In the early research, Abdel-Rahman used different probes of the hot wire anemometer (using the Dantec 55D90) to calibrate the speed and yaw angle. The results show that the yaw parameters change with the variation of the yaw angle and flow velocity [4]. Pesmajoglou and Graham applied numerical simulation by a free vortex lattice model and experimental tests to sequentially analyze the time-varying characteristics of the aerodynamic force and momentum coefficient of the impeller in stable yaw, sinusoidal oscillation yaw, and random yaw by an incident turbulent wind, in which there are three wind turbines including WTMR (LS(1)-0412), Marlec (NACA 4415), and Howden(LS(1)-0412). This work had obtained some correlation laws between the turbulent flow and the blade vibration force [5]. Schreck and Robinson et al., pointed out that the stall angle of attack of wind turbine (S809) during yaw is somewhat shifted from that of static ones, in which the surface pressures along the wind impeller blade were obtained at multiple span locations utilizing the NREL Phase IV unsteady fluid measurement. The rotor turns at a constant 71.6 r/min and has a maximum rated power of 19.8 kW. The stall vortex has deformation and amplification effects, which leads to increased blade normal force [6]. Grant and Mo et al., studied the wake vortex trajectory in scenes with various conditions of Marlec (NACA 4415) turbine yaw and blade azimuth based on laser sheet visualization (LSV) techniques and used the wake vortex model, in which each blade of impeller had been divided into a number of span-wise elements especially with a finer element distribution for tip part, to predict the flow field. The results indicate that the experimental data and the calculated values are consistent [7].

With the deepening of the research on the yaw wake of the horizontal-axis wind turbine, Tongchitpakdee and Benjanirat et al., conducted experiments and numerical analysis on the aerodynamic performance of such National Renewable Energy Laboratory (NREL) Phase VI horizontal-axis wind turbine at four wind speeds and four yaw angles, in which an unsteady compressible 3D viscous flow model for NREL Phase VI rotor are solved using an implicit time-marching scheme, and a third-order accurate upwind scheme is used to compute the inviscid fluxes at the finite element faces. They found that the flow is fully attached, which agrees well with measurements using a Baldwin–Lomax model at wind speeds less than 7 m/s. When the flow is over 20 m/s, it is massively separated, and the related quantities become insensitive to turbulence model effects. Between 7 and 15 m/s, the situation is partially separated, which needs more sophisticated turbulence models and transition models [8]. Maeda et al., conducted experimental measurements on the surface pressure distribution of horizontal-axis wind turbine impellers with rotor diameters of 10 m and 2.4 m under different yaw angles, in which spanwise airfoil combinations are formed with DU91-W2-250, DU93-W-210, NACA63-618, and NACA63-215 from the blade root to tip. They pointed out that as the yaw angle increases, the power coefficient of the rotor decreases, and the relative normal force coefficient is also reduced, which is compared to the non-yawed condition [9,10]. Haans and Kuik et al., used experimental measurements to study the aerodynamic performance of the impeller (NACA0012, with two blades, rotor tip radius of 0.6 m and root radius of 0.18 m) at an open jet wind tunnel in Delft University of Technology and the corresponding wake structure in the yaw motion of a wind turbine. This study combined the improved blade element momentum (BEM) theory to obtain the correlation law between yaw misalignment and turbulent eddies [11]. Based on the research of [11], Sant and Kuik et al., analyzed the characteristics of the adhering flow and boundary layer separation flow on the blade surface of the NREL Phase VI wind impeller for the periodic change of the angle of attack during the yawed conditions and proposed a reconstruction method for the time-dependent angle of attack and induced velocity distributions by developing a free wake vortex model, which does not directly apply the airfoil data but the spatially distributed vorticity generated by the blades and used the measured local blade loads [12]. Jiménez and Crespo et al., applied the large eddy simulation (LES) model with periodic boundary conditions implemented in a CFD solver to study the wind turbine wake deflection in yaw, in which the calculated results are consistent with the experimental data and illustrated that the Reynolds number does not have any significant impact on this issue, as the roughness of blade surface is the dominant factor in the flow, and proposed a new possibility for active wake controlling by using such simple analytical correlations in the study [13].

Following the development of experimental equipment and the advancement of computer technology, a series of new enhancements have been made in the investigation of the wake of wind turbines in yaw. Narayana and Putrus et al., designed a general simplified model to study the yaw behavior of small horizontal-axis wind turbines (NACA4415, wind rotor radius of 1.105 m, two blades) with a tilt-up impeller using a developed mathematical model based on D’Alembert’s principle, which indicates a dynamic system can be converted into a static system to be analyzed under static equilibrium conditions. They pointed out that when the wind speed exceeds the rated value, the wind rotor will tilt upwards and is facing the wind direction. The inclined position is stabilized, thereby controlling the wind component responsible for power generation. In this model, the aerodynamic effect is evaluated by the blade element theory (BET), and the induced velocity in the flow field is specifically predicted by the combination of the modified BET and the momentum theories [14]. Micallef and Bussel et al., used experimental and numerical methods to investigate the near-wake flow of a horizontal-axis wind turbine (MEXICO rotor, and root of DU91-W2-250, midboard of Risø A1-21, tip region of NACA 64418) using stereo particle image velocimetry (SPIV) and a 3D unsteady potential-flow panel model with the working conditions of axial and yawed flow (yaw angle from −30° to +30°). The results showed that the complex flow area near the root and tip in the axial state is related to blade vibration, and the radial velocity is not the main factor causing the vorticity change. In the yaw state, the tip vortex moves away from the blade quickly in the windward region, and the direction of the tip vortex movement in the leeward region is exactly the opposite [15]. Jeong and Kim et al., analyzed the effect of yaw error on the performance and dynamic stability of horizontal-axis wind turbine blades based on a free-vortex wake model and the corresponding measurements with the NREL Phase VI rotor of 5 MW and pointed out that yaw misalignment adversely affected the aeroelastic stability of horizontal-axis wind turbine blades [16]. Watanabe and Takahashi et al., derived a simplified equation of yaw rate which includes all design variables for the small wind turbine (NACA0012) and verified the validity of the equation through wind tunnel test results [17]. Afterward, Abedi and Davidson et al., developed a free vortex wake method coupling the theory of the potential, inviscid and irrotational flow to investigate the change of thrust vector using a two-bladed Hönö wind turbine [18]. Qiu and Wang et al., proposed an unsteady numerical simulation method based on a nonlinear lifting line method and a time-accurate free-vortex method to study the unsteady aerodynamic load of the blade and wake flow during yawing and/or pitching processes at different wind speeds for three different types of wind turbines (NREL Phase VI, TU Delft, Tjæreborg). It is pointed out that the blade loads, rotor torque, and the locations of the tip-vortex cores are well consistent with the experimental data, indicating that the novel method can give accurate aerodynamic loads with a better computational efficiency [19]. Bastankhah and Port’e-Agel used high-resolution particle image velocimetry (PIV) to measure the wake flow field of a wind turbine under different yaw angles and tip speed ratios and found that tip speed ratio and wake rotation are important factors that cause the highly complex wake structure in the near-wake region. The thrust and power coefficients of the impeller will decrease with the increase of the yaw angle, and the wake loss will increase with the increase of the yaw angle, and the mean streamwise velocity in the far wake in the condition of non-yaw angle displayed a self-similar Gaussian distribution [20]. Kress and Chokani et al., tested the yaw stability and performance of three different downwind rotors (with cone angles of 0°, 5°, and 10° based on a scale-model of the multi-megawatt impeller) through experiments and compared them with the corresponding upwind rotors. The results show that wind turbines with upwind have significantly lower yaw stability than wind turbines with the downwind rotor configurations. Under the same yaw angle, cone angles of 0°, 5°, and 10° downwind have higher shaft power and rotor thrust than the corresponding upwind configurations [21].

In recent research, Dai and Yang et al., gave the definition of yaw coefficient and its specific calculation method with SCADA data and subdivided the loss factor of energy capture (power coefficient) into aerodynamic loss factor and inertial loss factor in order to consider aerodynamics and mechanical inertia on energy capture. Studies have pointed out that if the wind speed remains constant with the rated value but the yaw coefficient increases, the corresponding output power will decrease. When inflow wind speed is close to the rated value or above it, the influence of the pitch control on the output power is stronger than the yaw [22]. Ke and Yu et al., studied the aerodynamic performance and wind-induced effect of a 5 MW wind turbine tower-blade system (developed by Nanjing University of Aeronautics and Astronautics with airfoils of the NH02 series) under different yaw angles (0°, 5°, 10°, 20°, 30°, and 45°) with the wind–rain combination action, revealing the velocity flow line, turbulence energy strength, raindrop running speed, and trajectory action mechanism in the coupling field, and constructed the large-scale wind turbine tower-blade coupling model with different yaw angles. Among them, the discrete phase model (DPM) in CFD simulation was added to achieve the wind–rain coupling synchronous iteration process. The results show that the yaw angle will significantly affect the aerodynamic and comprehensive stress performance of the wind turbine, and the wind–rain load enhances the structural response and reduces the overall buckling stability and ultimate bearing capacity of the system [23].

Saenz-Aguirre and Zulueta et al., proposed a yaw control strategy based on reinforcement learning (RL), considering multivariable states and actions, and analyzed the aerodynamic performance of the NREL 5 MW wind turbine with the aeroelastic code FAST. Among them, a particle swarm optimization and Pareto optimal front-based algorithm were improved to capture the optimal compromise-point between the power generation and the mechanical loads. They used real wind data in Salt Lake City, Utah, to verify the effectiveness of the designed yaw control strategy [24]. Saleem and Kim calculated and analyzed the effects of changes in yaw angle and wind speed on the performance of three horizontal-axis airborne wind turbines (NACA-5415, NACA-9415, and NACA-5425) with different shell shapes under high-altitude operating conditions, in which the k–ω SST turbulence model, based on Reynolds-averaged Navier–Stokes equations, was used. It is pointed out that the increase of the tip speed ratio and the yaw angle will increase the impeller thrust coefficient. In the large yaw angle condition, the instability of the buoyant shell increases. The shell based on NACA-9415 airfoil shows a higher power coefficient and stability [25]. Ye and Wang et al., used the sliding grid method and unsteady Reynolds-averaged Navier–Stokes computations to study the non-steady-state changes of the aerodynamic performance of small wind turbine (NREL Phase VI rotor with the airfoil of S809) during static yaw and yaw, and their 3D effect on the unsteady characteristics. Studies have shown that CFD results have better agreement with the experimental data than blade element momentum (BEM) results. The power of the wind turbine decreases as the yaw angle increases. The torque in the yaw state shows lower frequency fluctuations than the non-yaw state, and the fluctuation caused by the dynamic yaw is greater than the static yaw. At the same time, for dynamic yaw conditions, aerodynamic fluctuations are more pronounced on the back side than on the front side. As the yaw angle increases, the maximum load position moves from the midspan to the outside [26]. Kleusberg and Schlatter et al., employed the actuator-line method in the spectral-element code Nek5000 to conduct a numerical study on the wake development of a single yaw wind turbine operating with different tip speed ratios and yaw angles. The results show that, depending on the tip speed ratio, the blade load varies along the azimuth angle, resulting in asymmetric wakes in both the horizontal and vertical directions. A large tip speed ratio and a large yaw angle can reduce the vertical asymmetry of the counter-rotating vortex pair caused by yaw, and they increase the spanwise force caused by yaw. Particularly, the strength of the counter-rotating vortex pair in the far wake increases with the yaw angle, and the vertical change of the wake center highly depends on the yaw angle and the tip speed ratio [27].

In summary, for early and current research, including experimental measurements, numerical calculations, and theoretical analysis, it is shown that the changes in a series of wake flow characteristics caused by wind turbine yaw motion are mainly concentrated in static yaw angle conditions, and especially for the studies of wake vortices, parameters including the morphological change characteristics of the tip vortex and the central vortex and the corresponding quantitative are still unclear. More importantly, there are almost no reports on the three-dimensional morphological change characteristics and laws of wind turbine wake vortex system during the dynamic yawing.

Therefore, in the present paper, an unsteady computational fluid dynamics simulation against the continuous dynamic yawing of wind impeller based on the multi-relaxation time lattice Boltzmann method (MRT-LBM) is performed to investigate the dynamic changes of the 3D wake vortex system of a full wind turbine model. Additionally, considering the complexity of the flow during the dynamic yaw of the horizontal-axis wind turbine, the large eddy simulation (LES) model and the wall-adapting local eddy-viscosity (WALE) model are also used in this study.

The paper is organized as follows: Section 2 briefly describes the LBM methodology; Section 3 briefly gives a description of the large eddy simulation model; Section 4 describes the geometry model of the wind turbine, the spatial discrete layout of the square-lattice, and domain solving conditions; Section 5 provides results for the dynamic yaw wake characteristics and corresponding laws of change, as well as verification of the validity of the CFD model; Finally, Section 5 summarizes the conclusions drawn in this work.

2. A Brief Description of the LBM-LES Methodology

2.1. LBM Model

The LBM is a new CFD method that has emerged in recent decades [28,29,30,31], in which one of the particularly simple linear collision operators is proposed separately by Koelman [32] and Chen [33] using the Bhatnagar–Gross–Krook (BGK) collision model.

The lattice BGK model (LBGK) uses the local equilibrium distribution functions to satisfy the macroscopic fluid Navier–Stokes equations, but this type of model has a single relaxation time property, which shows irrationality in actual fluid dynamics. Afterward, an improved lattice Boltzmann method—namely, multi-relaxation time LBM (MRT-LBM), is proposed, which revised the definition of Prandtl number assignment and the ratio between kinematic viscosity and bulk viscosity [34,35].

The currently widely used lattice model was proposed by Qian and D’Humieres et al., in which the model is called the DdQm model with d-dimensional space and m discrete velocities [36]. Several typical lattice models in the two-dimensional and the three-dimensional space are shown in Figure 1.

In this study, the speed model D3Q19 was selected, and its model is shown in Figure 2.

This MRT-LBM model, which is more accurate than the LBGK model, is selected here, and its evolution equation is expressed as follows:

$$|f_{\alpha }\left({\mathit{r}}_i+{\mathit{e}}_{\alpha }\delta t,t+\delta t\right)⟩-|f_{\alpha }\left({\mathit{r}}_i,t\right)⟩=-\mathit{S}\left[|f_{\alpha }\left({\mathit{r}}_i,t\right)⟩-|f_{\alpha }^{eq}\left({\mathit{r}}_i,t\right)⟩\right]$$

(1)

where $f_{\alpha }$ represents the velocity distribution function, which satisfies the number density distribution containing kinetic particles at locations, $r_i$, and t denotes time. $e_{\alpha }$ is the discrete velocity vector, and $f_{\alpha }^{eq}$ is the equilibrium distribution function, and S is the collision matrix, which needs to be mapped through the corresponding conversion matrix for the actual calculation [35].

After adopting the D3Q19 model proposed by d’Humiéres, there are 19 discrete velocities on each lattice, which can be expressed as follows:

$$\begin{array}{cccccc}& & & & & \end{array}{e}_{\alpha }=\{\begin{array}{ccc}\left(\pm 1,0,0\right)\left(0,\pm 1,0\right)\left(0,0,\pm 1\right)& fo{r}_{}{}_{} \alpha =1,2,\cdots ,6;& \mathrm{group}\mathrm{A}\\ \left(\pm 1,\pm 1,0\right)\left(0,\pm 1,\pm 1\right)\left(\pm 1,0,\pm 1\right)& fo{r}_{}{}_{} \alpha =7,8,\cdots ,18;& \mathrm{group}\mathrm{B}\\ \left(0,0,0\right)& fo{r}_{}{}_{} \alpha =19;& \mathrm{rest}\mathrm{particle}\end{array}$$

(2)

The corresponding weight coefficient is

$$w_{\alpha }=\{\begin{array}{ccc}\frac{1}{18}& fo{r}_{}{}_{} \alpha =1,2,\cdots ,6;& \mathrm{group} \mathrm{A}\\ \frac{1}{36}& fo{r}_{}{}_{} \alpha =7,8,\cdots ,18;& \mathrm{group} \mathrm{B}\\ \frac{1}{3}& fo{r}_{}{}_{} \alpha =19;& \mathrm{rest} \mathrm{particle}\end{array}$$

(3)

Obviously, while ensuring good calculation accuracy, the D3Q19 model greatly reduces data storage, reduces the amount of calculation, and improves the calculation efficiency of numerical simulation.

Additionally, the scattering operator in the study is implemented in the central moment space, which can naturally improve the Galileo invariance, accuracy, and stability of the CFD code.

2.2. LBM-LES Hybrid Simulation Model

The hybrid calculation method is mainly to add the idea of LES to the computation process of LBM simulation. Dong and Sagaut pointed out that the overall equivalent viscosity is the sum of fluid molecular viscosity and turbulent viscosity [37].

$${\nu }_{}={}_{}{}^{}\nu {{}_0}_{}+{}_{}{}^{}\nu {}_t$$

(4)

where v₀ represents the fluid molecular viscosity coefficient, and v_t is the turbulent viscosity coefficient.

Since the relaxation time τ is defined as

$${\tau }_{}{=}_{}{\frac{\nu }{T}}_{}+\frac{\delta t}{2}$$

(5)

where T is the flow field temperature, and δ_t is the time step for unsteady computation.

Therefore, τ will become a changing function value instead of a fixed constant. In (4), the fluid molecular viscosity coefficient is given by

$${\nu }_0{}_{}{=}_{}\frac{UL}{Re}$$

(6)

where U stands for the velocity feature scale, and L is the length feature scale.

The turbulent viscosity coefficient v_t is defined by a sub-grid turbulence model of the LES model—namely, the wall-adapting local eddy (WALE) viscosity model, which provides a consistent local eddy viscosity and near-wall behavior [38].

$${\nu }_t{}_{}{=}_{}{\Delta }_f^2\frac{{\left(G_{ij}^d{G}_{ij}^d\right)}^{\frac{3}{2}}}{{\left(S_{ij}{S}_{ij}\right)}^{\frac{5}{2}}+{\left(G_{ij}^d{G}_{ij}^d\right)}^{\frac{5}{4}}}$$

(7)

$$S_{ij}=\frac{g_{ij}+g_{ji}}{2}$$

(8)

$$G_{ij}^d{}_{}{=}_{}\frac{1}{2}\left(g_{ij}^2+g_{ji}^2\right)-\frac{1}{3}{\delta }_{ij}{g}_{kk}^2$$

(9)

$$g_{ij}{=}_{}{\displaystyle \sum }_{\alpha =1}^{19}{e}_{\alpha i}{e}_{\alpha j}\left(f_{\alpha }-f_{\alpha }^{eq}\right)$$

(10)

where ${\Delta }_f=C_{\omega }\Delta x$, $\Delta x$ is the filter scale, and $C_{\omega }$ is the model coefficient, and the strain rate tensor g_ij is locally available with the non-equilibrium term of the distribution function of LBM. Here, the Cartesian lattice structure is well suited for the LES turbulence model due to the isotropic character of the turbulence out of the boundary layer.

For the important models of the sub-grid scale of LES and the processing model of the wall functions, readers can refer to studies in the literature [39,40,41,42,43,44].

3. Numerical Calculation Setup

3.1. Horizontal-Axis Wind Turbine 3D Model

The basic geometric model for CFD research is manifested in Figure 3, in which the sweeping surface diameter of the impeller is 1.4 m and the design wind speed of the wind turbine is 10 m/s. The rated power of a wind turbine is 300 W, and the starting wind speed is designed as 3 m/s.

In the process of 3D solid modeling of the horizontal-axis wind turbine impeller, we first adopted the reverse design method to extract and optimize the airfoil aerodynamic data of the sparrowhawk wing section. Then, it is standardized to form the standardized airfoil data format corresponding to the wind turbine blade. Finally, based on this kind of bionic airfoil data, the geometric model of the wind turbine impeller is established. The specific design parameters are shown in

Table 1. Operating parameters for the design wind turbine impeller.

Operating Parameter	Specific Value
Number of blades	3
Blade length	0.7 (m)
Diameter of impeller	1.4 (m)
Tip chord length	0.04 (m)
Ratio of span and chord	4.22
Tip twist angle	5.8 (deg)
Wind turbine airfoil	Sparrowhawk bionic airfoil
Design wind speed	10 (m/s)
Operating rated power	300 (W)
Starting wind speed	3 (m/s)
Tip modification	None

.

3.2. Computational Area Setting and Lattice Layout

According to the actual size characteristics of the experimental wind tunnel used by the cooperative team, we set a similar flow field computational domain, as shown in Figure 4, in order to better approximate the real flow situation of the wind turbine wake.

The geometry of the overall wind turbine has been arranged and placed centrally in the entire wind turbine domain, which occupies 9 (m) × 5.05 (m) × 7 (m) in the X, Y, and Z directions, respectively. Then, in order to we can better capture the detailed appearances of the tip flow area of the impeller, a cylindrical lattice encryption area is used to refine the grid discrete near the tip locations and the corresponding wake flow field area. The length of the cylindrical area along the X axis is 6 m, and its radius is 1.1 m.

The lattice space discrete structure is associated with the D3Q27 model, which adopts the spatial discrete form of an octree structure. It is a 3D scheme with the maximum number of discrete velocities in one computing element. The octree structure provides the possibility of using a non-uniform lattice structure, and therefore, different locations in the fluid domain have different spatial scales, in which the different spatial scales used are arranged hierarchically. The spatial–temporal scales of each level are twice smaller than the previous one, thus forming the above-mentioned octree structure. This is particularly effective because the dx/dt ratio remains constant throughout the fluid domain, thereby ensuring such constant CFL conditions and sound velocity at any location in the fluid computational domain.

For the resolution of the computational domain, considering the influence of spatial discrete lattice independence and the cost of computing resources, the global resolved scale is specified as 0.05 m, and the matching adaptive refinement algorithm is adopted to better capture the characteristics of the wake flow field of the dynamic yawing wind turbine. This high-precision computer-aided design (CAD) basic geometry file format is used, which can better characterize the complex aerodynamic shape of the wind turbine. In this way, the dynamic wake flow field of the wind turbine is more accurately approximated in the CFD calculation process. The target resolved scale for the reinforced cylindrical area is set to 0.025 m. In order to speed up the calculation efficiency of the fluid dynamics problem and enhance the stability of the numerical iteration in the calculation process, the discrete target scales of the non-encrypted area in the model is set to be consistent with the value of the global resolved scale, including series boundaries of the inlet, sidewalls, floor, and outlet. At the same time, considering the cell scale matching with the discrete area of the dense cylindrical space and improving the calculation accuracy of the localized flow area corresponding to the compound motion boundary of the wind turbine, we set that the surface target resolved scales of the impeller, tower, and tail rudder are all to 0.0125 m.

Consequently, the discrete layout of the computational domain is shown in Figure 5. The number of transition lattice levels is 3, in which there are different numbers of fluid cells in each sub-domain level. In the initial state of the computational region, the quantities of corresponding fluid elements are 1,856,268 in Level 1, 1,148,016 in Level 2, and 11,701,104 in Level 3, respectively. Due to the influence of the rotation and dynamic yawing motion of the impeller involved in the subsequent simulation, the number of transient computational elements in each sub-domain will change, but they are all controlled within a relatively small range. In a rigorous sense, the total number of fluid cells in the CFD model is tightly centered around the magnitude of 14.7 million.

The flexibility of the octree lattice structure and the advanced near-wall treatment applied by this study allows addressing one of the most important issues—namely, the fluid–structure interaction. We here handle dynamic impeller and overall yaw motion structure at the rear, involving the alternator and tail rudder, through the enforced behavior which moves the components based on one of the input positions and orientation laws and enforces the compound motion of the wind turbine.

For the above dynamic behavior, the lattice structure is updated every time step to mark the lattice nodes that belong to the fluid region and those that belong to the solid region. Additionally, the discrete velocities in each cell are also projected every time step in order to compute the new distance to the wall required for the wall function. The specific situation is depicted in Figure 6.

3.3. Determining the Solution Conditions

The inlet boundary of the CFD model is designated as velocity–inflow type, associated with the normal direction value of 10 m/s. The initial absolute pressure in the domain is considered a standard atmospheric pressure with the value of 101,325 Pa, and the level of initial turbulence intensity is assumed to be 0.5%. The type of outlet boundary is specified as the pressure-outlet boundary, in which the boundary pressure is consistent with the initial absolute pressure of the domain. The reference density of air is 1.225 kg/m³, and the corresponding operating temperature is set to 290 K.

The value of the design tip speed ratio of the wind turbine is 5.5, which corresponds to the impeller rotation speed of about 750.30 r/min, relative to the rotational axis of the wind turbine alternator. The yawing motion of the wind turbine is defined by a sine trigonometric function, which is used to characterize the periodic yaw swing of the impeller.

$$\vartheta =A\cdot sin((2\pi /T)\cdot t)$$

(11)

where A represents the absolute value of the limit of the yaw angle, and T equals the overall simulation time for the specified working condition in the study.

For the simulation of the turbulence, the CFD calculation is solved by the direct numerical solution in the LES at the large scale, and the WALE model is used to approximate the eddy viscosity within the sub-lattice scale, which depicts a coherent and coordinated local eddy viscosity, as well as flow characteristics near the wall. Additionally, in the CFD research, we also applied the wall-modeled LES approach (WMLES), in which the isotropic lattice structure will indicate an unsuitable high number for the cells to resolve the fluid behavior on the boundary layer.

This issue can be addressed by using a generalized law for the wall boundary.

$$\frac{U}{u_c}=\frac{U_1+U_2}{u_c}=\frac{u_{\tau }}{u_c}\frac{U_1}{u_{\tau }}+\frac{u_p}{u_c}\frac{U_2}{u_p}=\frac{{\tau }_{\omega }}{\rho u_{\tau }^2}\frac{u_{\tau }}{u_c}{f}_1\left(y^{+}\frac{u_{\tau }}{u_c}\right)+\frac{\frac{dp}{dx}}{\left|\frac{d{p}_{\omega }}{dx}\right|}\frac{u_p}{u_c}{f}_2\left(y^{+}\frac{u_p}{u_c}\right)$$

(12)

and

$$y^{+}=\frac{u_{c}y}{\nu },\hspace{1em}\hspace{1em}{{}_{}}_{}{u}_c=u_{\tau }+u_p$$

(13)

and

$$u_{\tau }=\sqrt{\frac{\left|{\tau }_{\omega }\right|}{\rho }},\hspace{1em}\hspace{1em}{{}_{}}_{}{u}_p={\left(\frac{\nu }{\rho }\left|\frac{d{p}_{\omega }}{dx}\right|\right)}^{\frac{1}{3}}$$

(14)

where y is the normal distance from the wall, and x is the local flow direction tangentially to the wall. The $u_{\tau }$ is the skin friction velocity, ${\tau }_{\omega }$ is the turbulent wall shear stress, $\frac{d{p}_{\omega }}{dx}$ is the wall pressure gradient, $u_p$ is a characteristic velocity of the adverse wall pressure gradient, and U is the mean velocity at a given distance from the wall, and the f₁ and f₂ are the interpolating functions, which are submitted to the unified laws of the wall boundary.

The physical velocity field near the wall boundary is approached by the y⁺ that depends on the distance between the first lattice from the geometry wall y and the velocity u_c for the first lattice.

4. Numerical Results for the Wake in Dynamic Yawing Process

In the unsteady CFD numerical computation, we adopted a relatively small advancing time step, namely about 8.48764 × 10⁻⁵ s, to ensure the stability and accuracy of the initial flow field calculation, which is very important for capturing the subsequent wake feature changes. The value of A in (11) is set to 30 deg according to the conditions of the corresponding experiment, and T which stands for the overall simulation time is 9 s.

4.1. Overall Arrangement of Selected Snapshot for the Unsteady Flow Wake

To clearly demonstrate the relevance of the selected wake snapshots during the dynamic yaw process, we give the relationship between the deflection angle and time and indicate the position of the selected observation points in the overall evolutionary advancement. In Figure 7, the red curve with the form of a sine function represents the dynamic continuous change process of the deflection angle of the wind turbine in the time domain. The series of blue dots denote the selected wake snapshots, which are selected according to the change of approximately every 10 degrees in deflection angle. In addition to the overall change characteristics and selected snapshots of dynamic yaw, more detailed information on deflection angle and time is also displayed.

4.2. Changing Characteristics in Near Wake of Impeller for Unsteady Yawing Flows

In the progress of calculation advance, a series of transient snapshots of the wake flow field of wind turbine impeller are obtained through the corresponding transient vorticity iso-surface with the value of 150 s⁻¹, according to the arrangement illustrated in Figure 7.

To be more conducive to the analysis and comparison of the corresponding results, we normalized the relevant physical quantity, in which we used different normalized processing methods for velocity and static pressure, due to their respective data characteristics. The specific normalization formula is as follows:

$$\hat{X}=\frac{\left(X-Min\left(X\right)\right)}{\left(Max\left(X\right)-Min\left(X\right)\right)}, \widehat{X}\in \left[0{,}_{}1\right]$$

(15)

where X represents the velocity vector, and $\widehat{X}$ is the normalized result for the flow field.

$$\widehat{Y}=\frac{Y}{Max\left(abs\left(Y\right)\right)}, \widehat{Y}\in \left[-1{,}_{}1\right]$$

(16)

where Y represents the static pressure vector, and $\widehat{Y}$ is the normalized result in the domain.

The snapshots from 1 to 13 are manifested one by one in an appropriate rendering mode, shown in Figure 8, and the important details of corresponding advance time and deflection angle are also presented. The initially typical localized flow structure around the impeller is effectively captured, which is three spiral wakes that extend backward and are approximately symmetrical in the circumferential direction, according to the captured snapshots of Figure 8a,b. The tip vortex is essentially a local flow around the tip caused by the pressure difference on the tip surfaces between windward and leeward, which is the existence of relatively large-scale eddies periodically formed due to the dynamic static pressure changes on the rotating blade surface. The initially extended spiral wake in Figure 8b has relatively strong local random fluctuation characteristics because the flow is at a high Reynolds number level with approximately the order of magnitude of 1.5 × 10⁶, in addition to the effect of dynamic yaw motion.

As the initially formed tip vortex contains relatively high rotational kinetic energy corresponding to large-scale vortex structure, and the influence of inertial force dominance in the high Reynolds flow status is significant. The localized rotating flow characteristics need to go through a certain period of the dissipative process to attenuate to the isotropic small-scale vortex level, which causes a distinctly continuous tip spiral wake structure, as in Figure 8c,d.

In accordance with the change characteristics of the deflection angle in Figure 7, the absolute value of the rate of change of the yaw angle from snapshot 1 to 2 is at a relatively high level, which makes the tip wake very unstable in Figure 8b. It is probably due to the rapid change of the aerodynamic angle of attack at the tip of the blade and the increased unevenness of the turbulence in the near wake. As the rate of change decreases slightly from snapshot 2 to 3, the stability of tip spiral wake has been improved to a certain extent as in Figure 8c. In the process from the snapshot 3 to 4, the rate of change of the yaw angle decreases further until it reaches zero, and the corresponding stability of tip spiral wake is also further enhanced even at the deflection angle with the largest absolute value of 30 (deg) in Figure 8d. A similar feature also exists in the processes from snapshot 4 to 7 and from snapshot 7 to 10, as well as from snapshot 10 to 13. It illustrates the degree of stability of tip spiral wake in the dynamic yaw condition is inversely related to the absolute value of the change rate of yaw angular speed. However, for the relatively large changes of the incoming wind, the relevant characteristics of the wake flow will be more complicated, owing to the fact that the alters of wind speed and direction will lead to increased the local variability of the angle of attack in aerodynamics along the blade span.

To further study the characteristics and laws of the complex near-wake flow of the wind turbine, we set up a series of feature surfaces to investigate the flow structures of the tip wake region and the center wake area in detail, as shown in Figure 9. This series of feature surfaces include longitudinal feature surface of Z = 0 (m), and transverse feature surfaces of X = 0 (m), X = 0.3 (m), X = 0.6 (m) and X = 0.9 (m), in which the feature surface of Z = 0 (m) is one of the most typical and effective ways to study the wake characteristics of wind turbines.

Consequently, we give a series of numerical computation results of wind turbine wake flow on the feature surface of Z = 0 (m) for the dynamic yaw process in Figure 10, according to the planning of the selected wake snapshots as in Figure 7. In Figure 10a, the initial swirling flow wake of the impeller is displayed with a yaw angle of 0.4 (deg). As the continuation of the compound motion of rotation and yaw, at the time of about 0.49 (s) with the yaw angle is approximately 10 (deg), the typical wake flow sub-region has been preliminarily formed, consisting of the tip vortex region, central wake area, and interaction region of tip vortex and tower shadow wake flow.

When the yaw angle reaches about 20 (deg), corresponding to snapshot 3 in Figure 10c, the position and shape of the tip vortex are more distinct in the near wake, and the speed loss zone in the central wake area is basically formed, corresponding to the blue area in the center in Figure 10c. The radial outward movement of the tip vortex is captured, but the symmetry of the tip vortex in the upper and lower parts of the near wake is affected to a certain extent, mainly due to the interference of the tower shadow effect. As the wake propagates downstream, because of the mutual fusion between the different vortex structures, the degree of velocity difference in the flow field decreases for the areas inside of each sub-region and between different sub-regions.

To the yaw angle of 30 (deg) in Figure 10d, the radial outward movement of the blade tip vortex shows the linearization, and the speed loss zone extends downstream for the central wake flow. Subsequently, the yaw angle from 30 (deg) to 10 (deg) in Figure 10d,f, the linearization of tip vortex continues, and the speed loss zone continues to expand downstream with more serious losses. When the wind turbine returns to the equilibrium position, as in Figure 10g, with the yaw angle of 0 (deg) approximately, the linearized migration of tip vortex is changed, which corresponds to the maximum angular velocity change rate as in Figure 7 and also is consistent with Figure 8g. The speed loss of the center wake is reduced at about the normalized velocity of 0.27, and the influence distance of the loss zone downstream is shortened. At the same time, the speed loss zone in the center of the near wake has a transverse expansion. It shows that the high acceleration of yaw angular will have a great impact on near wake both for the regions of tip and center, which is most likely because the aerodynamic angle of attack changes too drastically along the blade span at the high acceleration level. From snapshot 8 to snapshot 13, it exhibits similar changing characteristics, corresponding to Figure 10h,m.

According to the above calculation results as in Figure 10, the main change characteristics of the wake of the wind turbine are shown in the near-wake area; therefore, the subsequent transverse feature surfaces are mainly distributed in the near-wake region. As in Figure 11, the wake flow results of the feature surface of X = 0 (m) are manifested from snapshot 1 to snapshot 13.

The approximately symmetrical distribution of the near-wake vortex structure along the circumference is observed both for the regions of tip and center in Figure 11a,g,m, corresponding to the equilibrium position of the dynamic yaw compound motion. In the feature surface of X = 0 (m), the directional induce of the normal direction of impeller sweep surface for tip vortex structure is clearly reflected in Figure 11a–g with the sine change of the yaw angle, in which there is obvious tip vortex distribution on the entering side, and few on the exiting side, especially for the deflection angles with the largest absolute value, such as Figure 11d,j. Similarly, the directional inducing downstream for the tip vortex is also captured in Figure 11g–m, in which their three-dimensional shapes are related to Figure 8g–m, respectively.

To further study the evolution of vortex structure at different X-axis positions at the same snapshot, the congruent rectangular window along X-axis is used in the transverse feature surfaces of X = 0 (m), X = 0.3 (m), X = 0.6 (m) and X = 0.9 (m), as shown in Figure 12. In Figure 12a, as a result of the insufficient advancement of computation time at snapshot 1, the propagation of the wake vortex has not reached the subsequent feature surfaces of X = 0.3 (m), X = 0.6 (m) and X = 0.9 (m). According to snapshot 2, the situation of the wake propagation can be observed in the corresponding follow-up feature surfaces in Figure 12b. The directional inducing by the impeller sweep surface for tip vortex is weakened with spreading downstream, as shown in Figure 12c–f. Specifically, as the X-axial distance between the impeller sweep surface and the corresponding transverse feature surface increases, the tip vortex distributions appear both on the entering side and exiting side. This phenomenon is probably due to the coupling effect by the main wind direction and mutual fusion of wake vortex structures.

Additionally, the severity of the speed loss is increasing from Figure 12c–f, and the scope to be affected in the central wake region is also expanding but to the situation with non-yaw; the speed loss has recovered to a certain extent in the central region, as in Figure 12g. It implies that the wind turbine impeller is in a state of relatively high energy conversion efficiency at snapshot 7 with the approximately non-yaw condition. From snapshot 8 to snapshot 13 in Figure 12h–m, the propagation downstream of the near wake flow field has a similarly regular pattern.

4.3. Static Pressure on the Blade Surfaces

Here, closely related to the developmental wakes shown in the specified instant photos in Figure 8, Figure 10, Figure 11 and Figure 12, the corresponding snapshots of static pressure distribution on the blade surfaces both for pressure and suction sides are manifested in Figure 13. The speed and aerodynamic load of the blade tip of the horizontal-axis wind turbine are at a relatively high level. When a blade with a finite length rotates, due to the pressure difference between the tip pressure side and the suction side, the airflow from the pressure surface bypasses the tip surface and flows into the suction surface, which destroys the local flow pattern of the tip, and at the same time, the localized flow behavior generates tip vortices. Especially for the section from mid-span to tip region, the variation of the distribution of static pressure values on the blade surface is sensitive and important, caused by the varying yaw angle.

Due to the influence of dynamic yaw motion, in the initial stage of impeller rotation, the pressure side and suction side of the blade show relatively uniform static pressure distribution, which means that the aerodynamic characteristics of the impeller surface are delayed, shown in Figure 13a. To snapshot 2, as in Figure 13b, the relative motion between the impeller blade and air makes the surface of the blade produce a local airfoil influence difference distribution for static pressure value, such as the tip part marked in the black ellipses.

From snapshot 2 to snapshot 4, the influence of changes in the aerodynamic angle of attack dominates along the blade span, resulting in that the static pressure distribution affected by the local airfoil at snapshot 2 has also changed, particularly marked in the red ellipses, as shown Figure 13c,d. When the yaw angle is reduced in Figure 13e,f, the absolute value of the maximum negative pressure on the blade surface increases, especially the parts from the mid-span to the tip.

As to snapshot 7 at the state of non-yaw condition, the effect of the local airfoil again occupies the main role along the span to form the differentiated static pressure distribution, such as marked in the black ellipse in Figure 13g. This also shows that the periodic dynamic stability of the static pressure on the surface of the impeller is increasing compared to snapshot 1.

From snapshot 7 to snapshot 12, corresponding to Figure 13g,l, in such coupling effect between the dynamic yaw and the local airfoil, the influence of the former is weakening, and the latter is strengthening. Especially, when it comes to snapshot 10 to snapshot 12, the typical characteristics of the leading edge and the suction surface basically appear, which are marked in the black ellipses. For snapshot 13 in Figure 13m, the static pressure distribution characteristics of the leading edge, pressure surface, and the suction surface of the blade are more obvious, and the periodic dynamic stability of the static pressure on the blade surface is further enhanced, contrasting with snapshot 1 and snapshot 7. However, due to the influence of yaw angular acceleration with the highest level at snapshot 13 displayed in Figure 7, the static pressure distribution has a local fluctuation from mid-span to tip in the black ellipse. This point indicates that the coupling effect has always existed in the dynamic yaw process, even for those snapshots with a non-yaw angle.

4.4. Overall Aerodynamic Performance of the Horizontal-Axis Wind Turbine Impeller

The numerically computational time-advancement histories of the axial force and torque along the X direction are clearly plotted in Figure 14a,b, respectively. To check the effectiveness and accurateness of the calculated CFD model for wake flow in the dynamic yaw process, the experimental equivalent measurement results of the corresponding force and torque against X direction are also manifested in Figure 14.

In the wake flow calculation of the CFD model, the overall air force and torque along X direction are integrated on the surface of the impeller by an interval of each time step, which relies on a series of lattices that are organized in an octree-structure algorithm. The specific boundary setting of the calculation model and related theory refers to the content of Section 3.

Due to the fact that there are more aerodynamic influence factors and hardware transmission errors in the experimental measurement process, the overall value of the experimental data is larger than computation results and the fluctuation range is also wider both for the force and torque along the X direction. However, the periodic change frequency and the overall evolution trend, between calculation results and experimental data, are both consistent with each other. Therefore, this implies that the CFD model of wind turbine with LBM-LES is effective and adaptive.

According to the results of Figure 14a, as flow time progresses, the fluctuation range of the axial force shows a certain degree of convergence both for calculated and experimental situations, which indicates that the flow stability of the impeller wake is increasing, even under the influence of dynamic yaw. Additionally, it is also worth noting that the quasi-periodic fluctuation of axial force has a discontinuous interval between 4 s to 5 s, which is probably due to the adaptive recovery of the overall aerodynamic angle of attack of the blade surface of the impeller near around the yaw balance position.

For the change of the overall aerodynamic torque of the impeller along the X direction in Figure 14b, its main change characteristic trend is similar to the impeller aerodynamic force, which also shows a certain degree of convergence for the value of aerodynamic torque. Moreover, the central value of the torque fluctuation presents a very obvious trigonometric function change characteristic, which means a very strong correlation with the corresponding dynamic yaw motion mode defined by formula (11) in Section 3. Therefore, it can be considered that the aerodynamic torque is more sensitive to the air force, and the change rule of the torque average value of the overall fluctuation will largely depend on the traits of impeller compound motion and the characteristics of yaw angular velocity, especially for the typical reciprocating periodic yaw motion situations.

5. Conclusions

The research on the complex wake of horizontal-axis wind turbines during the dynamic yaw process is currently scarce, and for this reason, we performed an unsteady computational fluid dynamics simulation using the MRT-LBM method to investigate the dynamic changes of the near wake of wind turbine impeller under a yawing case. Moreover, considering the complexity of this wake flows with dynamic yaw, the wall-adapting local eddy-viscosity model was also used in this research. The unsteady complicated flow situation in the near wake was studied in detail, including the changing wake characteristics, the variation of static pressure distribution, and the overall aerodynamic performance of the impeller.

The results indicate that the degree of stability of tip spiral wake in the dynamic yaw condition is inversely related to the absolute value of the change rate of yaw angular speed. However, the relevant characteristics of the wake flow will be more complicated for the relatively large changes of the incoming wind, because of the increased variability of the angle of attack in aerodynamics along the blade span. When the wind turbine returns to the equilibrium position with the yaw angle of 0 (deg) approximately, the linearized migration of the tip vortex is changed, and the speed loss of the center wake is reduced and has another transverse expansion. The directional inducing downstream of the normal direction of impeller sweep surface for tip vortex is clearly reflected, in which there is obvious tip vortex distribution on the entering side, and few on the exiting side, especially for the deflection angles with the largest absolute value. This phenomenon is weakened with wake spreading downstream, probably due to the coupling effect by the main wind direction and mutual fusion of wake vortex structures.

Additionally, the periodic dynamic stability of the static pressure on the blade surface is continuously enhanced, and the coupling effect between the dynamic yaw and the local airfoil has always existed in the dynamic yaw process, even for those snapshots with a non-yaw angle. The aerodynamic torque of the impeller is more sensitive to the air force, in which the torque average value will largely depend on the traits of impeller compound motion and the variation of yaw angular velocity.