Direct numerical simulations of aerodynamic performance of wind turbine aerofoil by considering the blades active vibrations

In the present study, the aerodynamic performance of the horizontal-axis wind turbine blades by considering the ﬂ ap-wise oscillations are numerically investigated by using direct numerical simulations (DNS). The details of ﬂ ow structure can be analysed and predicted by performing DNS over an oscillating blade by considering the realistic behaviour of the wind turbine blade structure with natural vibration frequencies. In this study, the impact of vibrations on the ﬂ ow separation point, laminar separation bubble (LSB) and stall over NACA-4412 aerofoil are investigated utilising the high- ﬁ delity spectral-hp element methodology. The Reynolds number and angle of attack were selected in the range of 50 ; 000 (cid:1) Re (cid:1) 75 ; 000 and 8 (cid:3) (cid:1) AoA (cid:1) 16 (cid:3) . It is found that the blade vibrations have a noticeable impact on the aerodynamic performance and delay the stall occurrence, and the lift remains high even at higher AoAs, in comparison with the stationary blade. The size of the ﬂ ow separation is reduced by the blade oscillation and the vibration also affects the separation point. Due to the harmonic oscillation of the blade, the pressure signals are periodic, and the pressure ﬂ uctuations are ampli ﬁ ed by the oscillations, especially in the ﬂ ow separation region. The time-averaged lift coef ﬁ cient is increased by 255.3% by raising the angle of attack, from 0 (cid:3) to 12 (cid:3) at Re ¼ 75,000. Compared to Re ¼ 50,000, the peak-to-peak amplitude for the angle of attack of 0 (cid:3) is higher, whereas that of 8 (cid:3) and 12 (cid:3) are slightly lower at Re ¼ 75,000. © 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Literature review
With the reduction of fossil fuels in the near future, and the growing interest in expanding economic affordable clean energy resources, wind turbines, as one of the leading renewable energy production equipment, has been widely utilized in the past few years. Wind turbines (WT) have various shapes, and each one of them has its own advantages. Horizontal Axis Wind Turbines (HAWTs) are one of the most popular wind turbines, and they have many benefits compared to the other types of wind turbines [1]. However, one of the main challenges in the design of HAWTs is to control the flow behaviour on the blade surface, which has a considerable effect on the performance of the wind turbines and in the noise generation and dynamic stall. Consequently, it is essential to perform precise numerical simulations to investigate the flow separation point and dynamic stall by considering the blade vibrations, especially at low Reynolds numbers ð10 4 < Re < 10 5 Þ, suitable for urban applications of HAWTs, which were not fully investigated.
Several experimental and computational aerodynamic analyses have been performed over horizontal axis wind turbine blades in the past decades. Wahidi et al. [2] experimentally examined the flow structures and laminar separation bubble (LSB) formation over NACA 4412 WT blade by using the PIV visualisation technique. They found that the location of flow detachment and reattachment depends on the value of AoA. Moreover, the visualisations revealed that the pairs of vorticities were created in the separated area, and they have an influence on the pressure fluctuations and flow disturbance in the reattachment region. Naung et al. [3] numerically investigated the impact of instabilities on separated region and the pressure coefficient over the wind turbine blade. The simulations showed that the inflow wake has a noticeable influence on the flow detachment and vorticity production over the blade. The authors [4,5] applied the frequency-domain method to calculate the force response and flutter instabilities on the HAWT blade. It was deduced that by utilising the proposed solver, the computational time would be decreased 200% in comparison with the other methods. The impacts of the internal blades phase angle on flow detachment and recirculation flow in the separated region were also examined. The stall phenomenon at low airspeeds can be categorised into three groups: trailing edge stall, leading-edge stall and thin aerofoil stall or LSB [6]. Gaster [7] investigated the formation of LSB over airfoils at different Reynolds numbers. It was concluded that the size of the bubbles depends on the pressure gradient and separation point over the surface of the airfoil. Koca et al. [8] experimentally examined the vortex production over NACA-4412 WT aerofoil. They evaluated the impact of Re number and AoA on recirculating vorticity generation. In addition, they employed some visualisation techniques to observe the flow path lines on the aerofoil surface. O'meara and Mueller [9] performed an experimental study on LSB generation on the surface of NACA 663-018 aerofoil at ð5 Â10 4 Re 2 Â10 5 Þ and ð8 < AoA < 12 Þ. It was concluded that the size of LSB was diminished by increasing the Reynolds number. Previous numerical simulations were mostly concentrated on using low-fidelity numerical methods (Such as RANS models) to investigate the flow separation, LSB, and stall occurrence over turbine blades. Zhang et al. [10] employed the Reynolds-averaged Navier-Stokes (RANS) model for the numerical analysis of the flow behaviour over WT blade at stall conditions. The simulations revealed that the Strouhal number remained in the range of 0.16e0.2 at various sections of the model. In the computational work of Lee et al. [11], the authors used the non-linear vortex lattice method (NVLM) to investigate the impact of wake geometry on HAWT yaw aerodynamic performance. They found that a large yaw misalignment causes non-sinusoidal oscillations of the turbine blades. Zhang et al. [12] numerically analysed the performance improvement of different wind turbine shapes to find the optimised geometry for vertical-axis WTs. The computations showed that NACA 0018 airfoil with a pitch angle of 6 provides the highest aerodynamic performance among the tested geometries. They observed that the flow detachment happens faster on aerofoils with smaller thickness compared to thicker ones with increasing the AoA. Alam et al. [13] numerically investigated the LSB and the beginning of the transitional structure on the WT blade. The simulations showed that the key source of the flow instabilities is a convective source in these separation regions. However, utilising high-fidelity unsteady models are required to precisely study the flow behaviour and instant vortex generation on the horizontalaxis wind turbines.
Accurate prediction of the lift and drag coefficient by considering realistic physical conditions of the turbine blades is essential to analyse and improve the performance enhancement of wind turbines. Therefore, in the past few years, wind energy scientists have had great interest to perform direct numerical simulations (DNS) or large-eddy simulations (LES) and other high-fidelity techniques [14]. It will help designers to enhance the performance of the wind turbines and delay the stall by considering all the details of instantaneous flow characteristics at a wide range of Reynolds numbers. Naung et al. [15] provided a DNS study to calculate the interfaces amongst the transient flow and the blade vibrations of a turbine. They deduced that the harmonic balance technique requires far smaller computational resources in comparison with the time-domain method. Cui et al. [16] numerically investigated the separated flow behaviour over WT aerofoil by using k-u-g turbulence integration method. It was deduced that the suggested method could provide more accurate results compared to LES with respect to experimental data. In the numerical study of D'Alessandro et al. [17], LES was employed to capture the flow separation point and stall occurrence over dimpled NACA 64014 WT airfoil. The simulations showed that the dimpled side can decrease the LSB volume if the dimples were located prior to the occurrence of LSB.
At low Reynolds numbers, boundary-layer separation and stall are very common over the aerofoil blades, which decrease the aerodynamic performance of the WTs. Considering the additional energy added to the flow by the turbine blades with vibrations can efficiently control the dynamic stall occurrence and increase the performance [18]. Several active and passive methods have been employed to control the flow over various aerofoil blades. Passive techniques such as synthetic jet actuators [19], and vortex generators [20] were employed in the past few years to control the stall over the wind turbine blades. Moreover, Meijer et al. [21] provide the generalised formulation of the local piston theory to account for the arbitrary motion of aerofoils. Acıkel and Genc [22] employed partially deformable membranes for controlling the flow detachment and separated shear layer on NACA 4412 WT blade at relatively lower Re numbers. The experiments showed that the LSB formation was diminished by utilising the membrane control equipment. Besides, the lift coefficient was improved noticeably. Saleem and Kim [23] used computational fluid dynamics (CFD) to evaluate the impact of the tip clearance to control the stall and separation point over the NACA 9415 aerofoil blade of a wind turbine. They observed that the blade tip gap could enhance the power parameter noticeably.
Considering the deformation of the blades due to vibrations is one of the popular active control methods [24]. With the active oscillations of the blade surface, energy is included in the near-wall fluid to remain more attached to the surface. However, it strongly depends on the Reynolds number and the angle of attack. Lei et al. [25] performed two-dimensional numerical simulations on active oscillations of E387 airfoil at the Reynolds number of 30,000. They mentioned that using DNS can more precisely calculate the boundary layer formation on the blade surface. However, due to the high computation costs, they only used the RANS model for their computations. They concluded that the laminar separation could be controlled over the various angle of attacks (AoAs), while the lift coefficient was increased noticeably. In the experimental study of Munday and Jacob [26], the authors considered active vibrations over the suction side of NACA-4415 aerofoil by changing the camber. They measured the location of the flow detachment for various Re numbers and AoAs by employing a visualisation technique. They observed that the separated region was noticeably decreased on the oscillating WT blade compared to the stationary one. Khatam et al. [27] numerically investigated the separation control over a NACA-4415 blade at low Reynolds numbers. They concluded that the actuation could significantly reduce the reparation area over the blade. The recent study of Jokar et al. [28] revealed that the flap-wise vibrations of the HAWT blades play a vital role in the design of wind turbines and further investigations must be performed to investigate the flow structure over vibrating WTs. Therefore, it is essential to perform DNS over threedimensional HAWT airfoil by considering the blades flap-wise oscillations. It will significantly help designers to accurately predict the lift and drag coefficient, aeroelastic behaviour and the flow separation point over the wind turbine blades under realistic physical conditions.
A summary of the previous numerical studies and experiments done on oscillating aerofoil blades are provided in Table 1. It is found that previous studies were mostly concentrated on experiments or two-dimensional numerical analysis, and there is no available DNS study on three-dimensional wind turbine blades by considering active oscillations.

Spectral-hp element method
The spectral-hp element method is a high-fidelity and accurate DNS method to investigate the instantaneous vortex generation, flow separation, and aerodynamic characteristics over different objects [42]. This method is much more computationally efficient compared to the other DNS methods. Moreover, this method, which is implemented in an open-source code named NEKTARþþ has novel abilities to consider the effects of vibrations and movement of the geometry on the flow structure and complex flow disturbance on the surface of the object. Moreover, it is possible to perform three-dimensional simulations over airfoil geometries using Fast-Fourier-Transform (FFT) technique, which significantly decrease the computation time. This method provides the ability to investigate both incompressible [43], and compressible flows at different inflow velocities [44]. Moxey et al. [45] defined a new convergency control parameter called spectral vanishing viscosity (SVV) to control the aerodynamic simulations over more complex 3D geometries under various boundary conditions and physical parameters with higher polynomial orders. The study of Bao et al. [46] showed that the vortex induced oscillations can alter the pressure distribution and recirculating flow over the flexible risers. Cassinelli et al. [47] used the spectral-hp element method to examine the flow separation and LSB formation over a stationary low-pressure T106 blade. They used FFT technique to perform DNS over the 2.5D (linear) T106 blade. The instantaneous variations of the boundary-layer profile, Q-criterion, turbulent kinetic energy, and flow separation point were thoroughly investigated in their numerical study. However, it is essential to perform a high-fidelity DNS study over three-dimensional wind turbine blade by considering the impact of the blades aeroelastic vibrations on the flow structure and separation point on the surface of the blades.

Objectives and novelty
According to the literature review and summary of the previous studies in the field of the aerodynamics of oscillating WT blades in Table 1, the knowledge gap can be summarised as follows: The first gap is related to a limited number of high-fidelity numerical studies over oscillating wind turbine blades. Nearly all the previous numerical studies were used the RANS model and other low-fidelity methods. However, performing DNS is essential to precisely calculate the details of flow structure and separation point by considering the blade oscillations. The second gap is linked to the lack of three-dimensional simulations over vibrating wind turbine aerofoils. Most of the previous papers performed two-dimensional simulations, and there is only one recent paper which performed 3D analysis on vibrating blades. They used the URANS method, which is unable to precisely predict the details of instantaneous flow unsteadiness and LSB over the blades. The third gap is relevant to the range of Reynolds number and the wind turbine aerofoil shape. NACA 4412 aerofoil is very popular in the design of wind turbines. However, there are no numerical or experimental simulations over oscillating wind turbine blades with this kind of geometry. The fourth gap is linked to the direction of the blade vibrations. Most of the previous studies were focused on pitching oscillations of the aerofoils. The recent modal analysis of Jokar et al. [28] revealed that flap-wise vibrations are very common in the HAWT blades. However, there is no numerical or experimental aerodynamic analysis over wind turbine blades by considering the flap-wise oscillations.
To address these gaps, the present numerical work concentrates on the following objectives: To perform high-fidelity direct numerical simulations (DNS) over wind turbine blades by considering the blades flat-wise vibrations.

Numerical study
was selected in the range of low-speed airflows ð50; 000 Re 75000Þ to investigate the LSB formation and separation point on the surface of the aerofoil. The angle of attack (AoA) is in the range of 0e16 . The details of the aerodynamic forces on the mid-section blade and the definition of the other parameters are provided in Fig. 2. For the aeroelastic simulations, the modal analysis should be performed prior to the DNS to find the natural frequencies of the blade to use in oscillating blade cases. To blade is made of composite material, approximated by the orthotropic material properties, as modern wind turbines are constructed by composite materials to decrease weight. The details of the physical properties are provided in Ref. [48].

Mathematical modelling
The unsteady and incompressible continuity and momentum equations in the vector form can be expressed as: where u, p and n are the velocity vectors (u,v,w), pressure, and kinematic viscosity. A mapping function is used to implement the oscillations to the surface of the wind turbine. The latest master version of the NEKTARþþ is used for the aerodynamic analysis on local HPC cluster using 64 CPU nodes. The spectral-hp element method decomposes the computation domain into subdomains as U ¼ ∪ e2ε U e . In this equation, U is the computational domain and U e is a sub domain of the model. The mapping function converts the local cartesian coordinate ðx 1 ; x 2 Þ into a reference coordinate system of ðx 1 ; x 2 Þ in the mapped space such that: More details of the domain decomposing method can be found in Ref. [49]. The cartesian coordinate can be used to define the moving body in the spatial domain as [50]: where 4 n is a 2D function which is defined as a tensor product of j p .
In the present work, the continuous Galerkin (CG) kernel is selected for the transient incompressible flow simulations on the horizontal-axis WT blade. Eq. (1) could be simplified as: where NðuÞ ¼ uðu:VÞ and L ðuÞ ¼ V 2 u are the non-linear and linear convection and diffusion terms of the Navier-Stokes equation, respectively. Based on the velocity-correction scheme, the advection term could be defined by: The Poisson and Helmholtz terms, can be expressed as: By time-integrating of the N-S equation and applying the advection, Poisson and Helmholtz terms, the momentum equation in nþ1 st time-step is discretized as [51]: In which J e and J i are the explicit and implicit terms. The above equation can be written as: form of the Poisson term is calculated as: In Eq. (11), V4 is the operating function on the computation domain. By computation of the unknown parameters, lastly, u can be calculated as [43]: A high-order outflow boundary condition is selected for present study. It improves the stability of the DNS model and can accurately capture the vortex generation in the wake region of the blade. Eq. (13) is utilized to calculate the high-order outflow pressure: where S 0 ðn:uÞ ¼ 1 is a step function, u 0 is the characteristic velocity, and d > 0 is a dimensionless constant value [52].
Eq. (14) was employed to evaluate the velocity vector as: Fig. 3 displays of the solution algorithm based on the spectral-hp element method and the convergency control of the advection, Poisson, and Helmholtz terms, respectively. A model analysis was firstly performed to find the natural vibration frequencies of a wind turbine blade. Afterwards, a mapping function is employed to impose the harmonic vibrations on the surface of the blade. Fig. 4 shows the details of the mesh used in the computational domain. Every element is divided into five sections (P ¼ 5) on each edge. Prism mesh is used near the airfoil surface to precisely capture the vorticity generation and flow separation point, and to resolve the laminar viscous sublayer near the no-slip walls. Ten layers of boundary-layer mesh with the minimum size of 0.004 m for the closest layer with a growing ratio of 1.2 is used to resolve the laminar viscous sub-layer effects, and to capture the flow separation point precisely. Tetrahedral mesh is employed for the remaining sections of the computational domain. A refined mesh is used in the wake region of the wind turbine airfoil to capture the  instantaneous vorticity generation in this region. To ensure the accuracy of the boundary layer mesh in spectral/hp element method, the following equations must be satisfied: x þ < 20, y þ < 1, and z þ < 10 [53]. The viscous length is L * ¼ n=u t and u t ¼ where t w is the wall shear stress on the surface of the wind turbine blade. x þ , y þ and z þ values remain in the acceptable range to precisely capture the flow disturbance near the walls. The highest values of x þ ¼ 7.21, y þ ¼ 0.82, and z þ ¼ 3.25 where captured close to the leading edge of the airfoil due to formation of boundary layer on the surface of the blade. Therefore, the wall þ values are in acceptable range to precisely detect the vortex production and separation point on the suction and pressure sides of the oscillating wind turbine mid-section airfoil. The CouranteFriedrichseLewy (CFL) number (CFL ¼ UDx=Dt) was measured during the computations at each time step to ensure the convergency of the results, the. To avoid divergency of the computations, the CFL number must be kept below one (CFL<1). However, based on experience, the spectral-hp element method requires CFL number to be less than 0.5 to avoid divergency during the computations. The analysis showed that Dt ¼ 2 Â 10 5 is appropriate enough to keep the CFL number less than 0.5 during the whole computations with the polynomial order of 9.
A grid independence study is performed prior to the main DNS simulations to find the most appropriate mesh for the simulations. The mesh sensitivity analysis is performed at Re ¼ 50,000 over the stationary NACA-4412 airfoil with AoA of 8 in Table 2. The grid independence study in the spectral-hp element method is usually performed on the polynomial order (P) to find the most appropriate mesh without significant increment in the computation times. The location of the separation point is one of the essential parameters in the aerodynamic analysis of wind turbines. Therefore, the X sep /C is compared for different polynomial orders to find the most appropriate mesh. The total number of elements in the domain based on 48 planes in the span direction is also provided. The results show that the polynomial order of nine has negligible deviation (0.04%) with P ¼ 11 in capturing the flow separation point. Therefore, P ¼ 9, which corresponds to the total number of 10,003,824 elements, is selected for further simulations.

Validation
Before performing further simulations, it is important to validate the accuracy of the numerical results with experimental data available in the literature. To validate the DNS results over stationary wind turbine blade at a specific Reynolds number, the drag and lift coefficient at different angles of attacks (AoAs) were compared with the experimental results of Koca et al. [8] in Fig. 5. It is observed that the C L and C D coefficients of the present study are in good agreement with the experiments at Re ¼ 50; 000. It is observed that the lift coefficient is increased by raising the angle of attack, and after AoA ¼ 14 stall occurs, and the lift coefficient suddenly decreased while the drag coefficient is increased.

Modal analysis
The modal coupling method is used in the present study for modelling and simulation of the fluid-structure interaction. Using this method, the mode shapes need to be imported into the fluid solver prior to the flow simulation. The mode shapes are used to define the blade vibration. To obtain the mode shapes of the WT blade, the modal analysis is performed using a Finite Element Analysis (FEA) method. The dimensionless displacement is presented as the ratio to the actual displacement to the maximum displacement value D Dmax . The first four mode shapes, along with natural frequencies, are presented in Fig. 6. The summary of the modal masses and effective translational masses for the first four modes obtained from the modal analysis are listed in Table 3. The fundamental mode, which is the first mode shape, is employed in this study as the vibration mode. The first natural frequency of the NACA-4412 wind turbine blade is 0.87Hz, which is relatively good compared to the DTU 10 MW wind turbine (0.64 Hz), is defined to be the oscillation frequency in this study. The blade vibration is defined using the following harmonic function based on the vibration frequency and amplitude. The amplitude of vibration is set to 1% of the chord length. Fig. 7 shows the behaviour of the blade vibration and the stress distribution due to the vibration. It is shown that the maximum stress concentration is found in the blade root and inner region, and it is distributed  along the blade span. In the blade mid-section, where the simulation is performed in this study, a strong stress concentration and distribution is still observed. Therefore, to ensure structural integrity, an intensive analysis needs to be performed to determine the impact of the blade vibration on unsteady flow structures.

Aerodynamic analysis
The effects of the blade vibration on the unsteady flow structures are discussed in this section. The vorticity structures generated over the wind turbine blade aerofoil at the angle of attack of 0 at Re ¼ 50,000 are depicted in Fig. 8. The vorticity contours from both stationary and oscillating aerofoils are provided to highlight the effect of vibration. When the blade is not vibrating, the flow remains attached to the pressure surface, whereas it separates on the suction side. The flow recirculation and separation bubbles are observed in the separation region on the suction surface. The separated shear layers roll up and then shed from the trailing edge (TE) of the blade. The vortex shedding from the TE is like the K arm an vortex street. However, in the case of the vibrating blade, it is seen that the size of the flow separation is slightly reduced. More intense separation bubbles are created in the separated region, and the vortex generation becomes stronger due to the blade vibration. It is noticed that the size of the vorticity structures left from the trailing edge is smaller than those of the stationary blade case. Therefore, these results indicate that blade vibration influences the flow detachment and reattachment process.
The instantaneous lift coefficients generated over the vibrating blade at three angles of attack of 0 , 8 , and 12 at Re ¼ 50,000 and 75,000 are illustrated in Fig. 9. This figure shows the variation of lift coefficients over time due to the blade vibration at different angles of attack and Reynolds numbers. It is seen that the lift is higher at greater Reynolds numbers. The fluctuation of lift coefficients over time is mainly related to the blade vibration, which also shows that the blade vibration has an effect on the generation of lift. At AoA ¼ 0 , it is clearly seen that the peak-to-peak fluctuation is stronger than other angles of attack, and the amplitude of peak-to-   peak variation of lift coefficient is approximately 0.4 at Re ¼ 50,000. This variation is diminished by expanding the AoA. The peak-topeak variations of lift coefficient for A0A ¼ 8 and 12 are 0.3 and 0.1, respectively. It is observed that this variation is reduced by around 55% when increasing the angle of attack from 0 to 12 . Similar behaviour is observed at Re ¼ 75,000. However, the peak-to-peak amplitudes are larger at this Reynolds number. At AoA ¼ 0 , the lift coefficient varies between approximately 0.04 and 0.56, and the average value is found to be 0.311, which is also greater than that of Re ¼ 50,000 (i.e., 0.2). The time-averaged lift coefficient at AoA ¼ 12 is 1.105. The lift coefficient is increased by 255.3% when raising the angle of attack, from 0 to 12 . The amplitudes of the peak-to-peak fluctuations for the angles of attack of 8 and 12 are 0.2 and 0.08, respectively. Compared to Re ¼ 50,000, the peak-to-peak amplitude for the angle of attack of 0 is higher, whereas that of 8 and 12 are slightly lower at Re ¼ 75,000. Overall, the results show that the blade vibration causes the variation of lift coefficient, and the fluctuations are stronger at lower AoAs at both Reynolds numbers. Consequently, it can be deduced that the stability of the blade is at risk at lower AoAs when the blade undergoes vibration. Raising the Reynolds number can increase the average lift coefficient. Fig. 10 demonstrates the velocity contours and streamlines around the oscillating aerofoil at 0 and 16 angles of attack at Re ¼ 75,000, which shows the effects of the blade vibration on the development of vortex structures at relatively low and high angles of attack. At both angles of attack, laminar separation bubbles (LSB) are formed within the separation zone on the suction surface of the aerofoil. The flow separation is much larger at the angle of attack of 16 as expected. Because of the larger BL separated region, the separation bubbles and flow recirculation developed in the separation zone are much stronger at this angle than in the 0 case. At both angles, the separated shear layer and flow recirculation rolled up and shed away from the TE of the aerofoil.
The blade vibration not only influences the vortex generation process but also affects the flow separation. Table 4 presents the flow separation points at different Reynolds numbers and angles of attack from the stationary and vibrating blade cases. In both cases, by increasing the angle of attack, the separation point moves towards the leading edge, whereas raising the Reynolds number delays the flow separation and moves the separation point towards the trailing edge. It is seen that the flow separation points are all affected by the blade vibration. At Re ¼ 50,000, the separation point shifts slightly towards the leading edge at the angles of attack up to 12 , while it shifts towards the trailing edge at angles greater than 12 . Similar behaviour is also seen at Re ¼ 75,000. However, at this Reynolds number, the separation point shifts towards the leading edge at the angles lower than 12 and it shifts towards the trailing edge at the angle of attack of 12 or greater. The results from the experiment are also added for the 8 angle of attack to indicate that the current simulations estimated the separation points accurately. Pressure signals within boundary layers at different axial locations of the aerofoil in terms of X/C are illustrated in Fig. 11. The results are presented for the angle of attack of 0 . To highlight the impact of blade vibration, pressure signals from both stationary and vibrating blade cases are compared. In the case of the stationary blade, the pressure signals fluctuate around the average value. These fluctuations can be considered as turbulent fluctuations, which can only be captured using high-fidelity simulations such as direct numerical simulations. At X/C ¼ 0.15, which is near the leading edge, relatively smooth pressure signals with some variations are seen, and it is related to the laminar boundary layer. At X/ C ¼ 0.4, which is it is close to the flow separation point, the start of laminar to turbulent transition can be identified as the pressure signals become stronger. By moving the axial location slightly towards the TE to X/C ¼ 0.55, it is observed that the fluctuation of pressure signals become much stronger. This is because of the flow detachment and the creation of pressure bubbles within the flow separation zone. On the other hand, in the case of the vibrating blade, the harmonic oscillations of pressure signals are observed due to the periodic vibration of the blade. At X/C ¼ 0.15, the pattern of pressure signals is smooth. As it moves towards the trailing edge and at X/C ¼ 0.4, which is in the vicinity of the flow separation, the turbulent fluctuations become present and noticeable. These turbulent fluctuations and disturbances are amplified at X/C ¼ 0.55 due to the formation of separation bubbles and vortex generation.   At a relatively high AoAs, (AoA ¼ 16 ), the flow detachment on the suction surface occurs much faster compared to the 0 case. As a result of the large flow separation, the separation bubbles and vortex generation become much stronger in the separation zone on the suction surface, which activates more disturbances and perturbation to the downstream wake region. It is clearly observed that the flow separation is dramatically reduced by the blade vibration at this high AoA. The flow recirculation and strong separation bubbles develop at the flow separating point and travel near the TE. While the shed is away from the trailing edge, the downstream vortex shedding is found to be highly turbulent. It can be concluded that the unsteady flow is highly distorted by the blade vibration, and the blade vibration primarily influences the vorticity structure and the location of the flow detachment and reattachment. Fig. 14 shows the comparisons of the lift and drag coefficients versus AoA at the Reynolds numbers of 50,000 and 75,000 for the stationary blade and vibrating blade cases. The time-averaged parameters were measured over 20,000 time steps. Technically, the lift coefficient rises as the angle of attack is increased before reaching its peak at the stall angle. At Re ¼ 50,000, the stall is identified at 12 when the blade is stationary. In the case of the vibrating blade, the lift is slightly reduced due to the vibration from 0 to 12 ; however, a significant difference between the two cases can be seen after 12 . The stall is mainly delayed by the blade vibration as the maximum lift is observed at 18 . Ever after the stall angle, the lift remains high compared to the stationary blade case. Similar behaviour is also observed at Re ¼ 75,000. Raising Reynolds number delays the stall in both cases. The stall angles are found at 14 and 19 in the non-oscillating airfoil and the oscillating blade, respectively. At AoA ¼ 20 which is in the post-stall region, the lift coefficient is higher than that of the stationary blade by around 40% at both Reynolds numbers. On the other hand, a similar trend of the drag coefficient is observed for non-oscillating and oscillating airfoils at both Reynolds numbers, but the drag is slightly higher in the oscillating blade. Hence, it is deduced that the blade vibration is crucial and has a significant impact on not only the unsteady flow structures but also the aerodynamic loads as the lift coefficient is strongly influenced by the blade oscillations.
Due to the blade vibration, the wind turbine blade experiences unsteady aerodynamic loads, resulting from unsteady pressure distributions. In an aeroelasticity analysis, unsteady pressure distributions are often expressed in terms of the unsteady pressure coefficient and phase angle, which are calculated based on the transient perturbations. Fig. 15 shows the variation of the pressure phase angle on the airfoil surfaces for the AoA of 8 and 16 . At AoA ¼ 8 , the pressure phase angle is distributed between 20 and 100 , whereas this variation is larger at AoA ¼ 16 , with the phase angle distributing between 20 and 16 . The latter is related to the large flow separation. Fig. 16 demonstrates the unsteady pressure coefficient variations on the airfoil at different angles of attack. At AoA ¼ 16 , the unsteady pressure difference between the two surfaces is greater near the leading edge due to the larger angle of attack, which continues until X/C ¼ 0.8. After this point, the pressure distribution is distorted, and the fluctuation becomes stronger due to the flow separation and strong vortex production. In contrast, the pressure drop among the blade upper and lower sides is smaller at a lower AoAs, and the fluctuations near the trailing edge are not as strong as that of 16 . Fig. 17 illustrates instantaneous pressure contours over both stationary and vibrating blades at different AoAs at Re ¼ 75,000, which indicates the impact of blade oscillations on the pressure contours. Generally, the pressure is low on the suction surface and high on the pressure surface. The pressure difference between the two surfaces generates the lift force needed for the rotation of the rotor. Due to the flow separation, pressure bubbles are developed on the suction surface of the aerofoil, and they shed away from the trailing edge. As the angle of attack rises and the flow separation gets larger, the evolution of pressure bubbles also becomes stronger. In the case of the stationary blade at AoA ¼ 16 , highly unsteady pressure fluctuations and a strong generation of pressure bubbles can be observed as a result of a relatively high angle of attack. The blade vibration further amplifies the pressure fluctuation and the development of pressure bubbles. At AoA ¼ 0 , it is seen that pressure bubbles are also developed on the pressure surface due to the rolling up of the vorticity structures, which was seen in Fig. 12.
At AoA ¼ 16 , it can be noticed that the flow separation is significantly reduced by the blade vibration, which is in agreement with the vorticity fields discussed above. Despite the smaller flow separation, a strong generation of pressure bubbles can be observed on the suction surface of the aerofoil, which indicates that the blade vibration highly influences the pressure distribution around the aerofoil and adds more disturbances to the flow by generating pressure bubbles. The effects of wind turbine blade oscillations on the wake profile at X/C ¼ 0.1 are presented in Fig. 18. The results are compared with the DNS results for the stationary blade and also with the experimental results of Koca et al. [8] at Re ¼ 75,000. In the experiments, the authors used hot-wire probes to capture the velocity distribution and the Reynolds stress in the wake region of the NACA-4412 airfoil at different inflow velocities and AoAs. The probes were located at 0.1C, and the 20480 sample data points at the rate of 2 kHz were collected in the experiments. The dimensionless velocity is dropped from 1.04 to 0.61 over the stationary blade. Besides, noticeable fluctuations are detected in the wake profile in the downstream region of the vibrating blade. The main physical reason for the fluctuations in the wake profile is related to the additional flow disturbance over the vibrating blade. Moreover, the DNS results for the stationary case have agreed well with the experimental results. It is observed that the wake peak value was increased by raising the angle of attack from 8 to 16 . As discussed earlier, the size of the separation region is much larger for the higher angle of attacks.
To examine the interaction between the fluid and the blade surface, the power spectral density (PSD) of instantaneous variations of pressure at a specific point near the trailing edge of the wind turbine aerofoil blade is presented in Fig. 19. The results are depicted for the oscillating blades with two different angles of attacks ða ¼ 8 and 16 Þ. fC/U is the Strouhal number, and f is the frequency. It is found that the numerical simulations for both cases are converged for the polynomial order of 11.

Conclusion
In this paper, high-resolution direct numerical simulations are conducted on an oscillating wind turbine blade aerofoil, and the effect of the blade vibration on the aerodynamic performances and the transitional flow structures are thoroughly investigated. The main conclusions drawn based on the results obtained are as follow: The unsteady flow and flow behaviour are highly distorted and affected by the blade vibrations. More intense separation bubbles are created in the separated region and the vortex generation becomes stronger due to the blade vibration. At AoA ¼ 0 , the lift coefficient varies between 0.04 and 0.56, and the average value is found to be 0.311, which is also greater than that of Re ¼ 50,000 (i.e., 0.2). The amplitudes of the peakto-peak fluctuations for the angles of attack of 8 and 12 are 0.2 and 0.08, respectively. The lift coefficient is increased by 51.7% when raising the angle of attack, from 0 to 12 . The size of the flow separation is reduced by the blade oscillation and the vibration also affects the separation point. The laminar separation bubbles are intensified by the blade vibration and as a result, the vortex generation becomes stronger. Lift and drag coefficients show that the stall is delayed by the blade oscillations, and the lift stays high even after the stall, compared to the stationary blade case. It is also found that the peak-to-peak variation of lift coefficients due to the blade vibration are much higher at lower angles of attack, which could affect the stability of the blade. Raising the Reynolds number can increase the average lift coefficient. Due to the harmonic oscillation of the blade, the pressure signals are periodic, and the pressure fluctuations are amplified by the oscillation, especially in the flow separation region. These unsteady perturbations are stronger at higher angles of attack.

Data availability statement
All research data supporting this publication are directly available within this publication.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.