Numerical Study of a Horizontal Wind Turbine under Yaw Conditions

Over recent years, considerable attention has been devoted to the optimization of energy production in wind farms, where yaw angles can play a significant role. In order to quantify and maximize such potential power, the simulation of wakes is vital. In the present study, an actuator linemodel code was implemented in the OpenFOAMflow solver. A tip treatment was applied to involve the tip effect induced by the pressure equalization from the suction and pressure sides.(e Leishman–Beddoes dynamic stall (LBDS) model modified by Sheng et al. was employed to consider the dynamic stall phenomenon. (e developed ALM-CFD solver was validated for the NREL Phase VI wind turbine reference case. (e solver was then used in simulating the yawed wind turbine, and power variation was compared with UBEM and CFD. Overall, according to the obtained data, the coupled solver compared well with CFD. (ere was an improvement in terms of prediction of the phase delay that is due to the dynamic stall. However, there was still negligible overestimation in deep stall conditions. Based on the obtained results, it is suggested that the reduction of power output follows a cosine to the power of X function of the yaw angle. In terms of visualizing wake, the results demonstrated that the current ALM code was satisfying enough to simulate skewed wake and vortices trajectory. (e effect of advancing and retreating blade was captured. It was found that yaw led to the concentration of the induced velocity downstream, resulting in a lower velocity deficit on a broader area, which is essential for wind farm optimization.


Introduction
e efforts towards reducing the environmental impact caused by fossil fuel-based electricity production led to the definition of more sustainable and green power generation systems, among which wind turbines play one of the most relevant roles. During the last decades, wind has gone from an alternative to a mainstream energy source, which is currently driving the renewables evolution. e amount of electricity production covered by wind points out relevant shares in many countries, providing 15% of EU demand in 2019 [1] and showing a world average of 4.8% in 2018 [2]. With an average growth of installed power at 13.4% in the last ten years and the levelized cost of energy (LCOE) dropping for both onshore and offshore applications [1,3,4], wind turbines gather strong research interest. However, unlike the vertical axis wind turbines (VAWT) that are wind direction independent [5,6], the effect of yaw angle in horizontal axis wind turbines (HAWT) is still an important issue.
In this context, the definition of numerical models to simulate wind turbines in different operating conditions is of utmost importance. Focusing on the aerodynamic analysis of the turbine wake, the actuator line method (ALM) is widely adopted, representing a valid compromise between low-order models and 3D calculations of the actual geometry. e method introduces a body-force scheme within a 2D or 3D flow solver. In ALM, the blades are replaced by rotating lines, along which a force distribution is imposed. is is determined through tabulated airfoil data. Since the turbine geometry is removed, no boundary layer resolution is required, yielding a lean grid and a much lower computational effort. Furthermore, if a wind farm installation is considered, the prediction of the wake interaction through an ALM implementation greatly simplifies the simulation complexity.
In recent years, several authors performed numerical analyses employing ALM. Mikkelsen et al. [7] employed an ALM combined with FLEX5 aeroelastic code and the LES solver EllipSys3D to study the wake interaction between three HAWTs in a row based on the Tjaereborg turbine geometry. Optimal pitch settings were investigated. e same flow solver was coupled to an ALM implementation by Shen et al. [8], who performed a validation considering Mexico and NREL Phase VI wind turbines under yaw conditions. Churchfield et al. [9] coupled an ALM code to an OpenFOAM LES solver. Improvements in velocity and body-force distribution function were introduced, and the code was validated on NREL Phase VI and NREL 5-MW wind turbines. Positive results were presented for the prediction of near wake behavior and tip vortices. Weihing et al. [10] developed an ALM running with the LES solver FLOWer. Two commercial wind turbines were considered under offshore and onshore conditions, and good agreement was achieved in comparison with 3D URANS simulations. Matiz-Chicacausa and Lopez [11] investigated the NREL Phase VI in downwind configuration, focusing on the tower shadow effect. An OpenFOAM URANS solver was employed, and the turbine was simulated using both ALM modeling and the full 3D geometry, pointing out matching with measurements. Wang et al. [12] studied the wake of two offshore NREL 5-MW wind turbines through ALM and OpenFOAM LES solver. Different yaw and tilt angles were tested, and a control strategy was proposed. Ravensbergen et al. [13] implemented an ALM within a variational multiscale flow solver, and NREL 5-MW, NTNU, and NREL Phase VI wind turbines were analyzed. e simulations were carried out in standalone, back-to-back, and complex terrain conditions.
Despite several benefits, the application of ALM in the literature has been limited to validation in yaw conditions [8]. In fact, further studies are required to reach more comprehensive and detailed understanding and clarify the weaknesses and powers of this method in yaw conditions. e present paper aims to investigate the application of ALM in the prediction of wind turbine wake in yaw operation, which is vital to achieve reliable wind farm design. e method is implemented within the OpenFOAM flow solver. A relatively new dynamic stall model that is modified by Sheng for wind turbine conditions [14] is employed. e model is validated on the NREL Phase VI test case. e results are then compared with a lower-order method and a higher-fidelity approach, namely, 1D BEM and 3D URANS. e outline of the remaining paper is as follows. In the next section, the computational model that includes the employed actuator line model, dynamic stall model, URANS governing equations, and geometric model are described. In the results and discussion section, after verification and validation of the current model, the result of yawed case is investigated. Finally, a summary and conclusions section is provided.

Actuator Line Model (ALM).
OpenFOAM package is a set of C++ libraries meant to solve partial differential equations. In the current paper, a C++ library was developed for the OpenFOAM toolbox to implement the actuator line model. e ALM employs a series of spanwise element points, representing sections with a constant airfoil, twist, chord, and oncoming wind. e volumetric forces caused by sections are projected into the flow field. is procedure is accomplished through source terms in the momentum equations and applies the same forces on the fluid as the turbine blades. e blade aerodynamic forces are calculated by coefficients from tabulated airfoil data. ese forces are formulated as follows: where α is the local angle of attack measured as the angle between the chord and the local relative flow, C l is the lift coefficient, and C d is the drag coefficient obtained from tabulated airfoil data. e lift and drag coefficients are linearly interpolated for the local Reynolds number. Here, U → rel denotes the relative velocity, resulting from inflow velocity (U → i ) and the element velocity (U → e ). A sampling approach was employed to interpolate U → i from the CFD resolved velocity field. A total of 16 sampling points were considered at a distance of 2Δx grid from the location of elements.
en, a Gaussian function is utilized for the smooth projection of the calculated forces onto the flow field. It is defined as follows [15]: where F field is the volumetric force field in cells' position, f tip is the tip correction function, and F → characterizes the forces at actuator element points that are calculated by lift and drag forces. In this equation, N is the body index (i.e., blades, hub, and tower), and e describes the actuator element point index. A tip treatment was applied to consider the tip effect induced by the pressure equalization from the suction and pressure sides [16].
where B is the number of blades, β is the pitch angle (including twist), R is the blade radius, r describes the location of the element point in the spanwise direction, and λ is the tip speed ratio.

Dynamic Stall Model.
In this study, the modified Leishman-Beddoes dynamic stall (LB-DS) model by Sheng et al. [14] is utilized. Since the LB-DS model was first proposed for helicopter applications [17], Sheng has modified the model to coincide with the use of HAWTs. In this situation, flow is typically incompressible, and airfoils are relatively thick. e principles of this model are presented in this section. For more details, readers are referred to references [14,18]. e LB-based models include three parts: the unsteady attached flow, stall onset, and separated flow. e air loads for the attached flow consist of impulsive and circulatory modules. e total normal force coefficient in the unsteady attached flow condition is equal to the sum of impulsive and circulatory normal force coefficient [14]: where H n is a deficiency function that is defined based on the LB 3G model [19], C N a is the normal force curve slope, and α E n is an equivalent of incidence, which is calculated by the geometrical angle of attack (α n ) and three deficiency functions.
e deficiency functions of X n , Y n , and Z n are empirically driven by using the flow velocity and pitch rate. Additionally, there is a lag due to leading-edge pressure and boundary layer development, which is implemented by applying a first-order lag to α n : where α n ′ is the delayed angle of attack, and the deficiency function D α n is introduced as follows [14]: where T α is a time constant associated with pressure response, and it is principally independent of the airfoil shape. e nondimensional time step is Sheng et al. defined a stall-onset criterion for the angle that the dynamic stall may begin. is criterion considers a reduced pitch rate [14]. α n ′ > α cr n ⟶ stall, α cr n � α ds0 , α ss + α ds0 − α ss r n r 0 , r n ≥ r 0 , r n < r 0 .
e reduced pitch rate is expressed by the following equation: e separated flow solution includes two indicial responses: (I) trailing-edge separation and (II) leading-edge vortex shedding. e trailing-edge separation module is associated with nonlinear effects of trailing-edge flow separation through a dimensionless separation point parameter.
is parameter is attained by Kirchhoff flow approximation in unsteady conditions. It is described as follows [14]: where α 1 , S 1 , and S 2 are constants that change with Reynolds number. In addition to the pressure response delay, boundary layer (viscous) lag is formulated as follows [14]: where D f n is the deficiency function, and it is defined as follows [14]: where T f is a time constant, which is dependent on the Mach number. e normal force coefficient for dynamic separation point is superimposed with impulsive loading as follows [14]: ere is an additional overshoot in the normal force during vortex shedding. e vortex-induced normal force coefficient is expressed as follows [14]: where f n ″ and f n are the dynamic and static separation points, respectively. V x is a modulation parameter, which is given by a periodic function that incorporates the track of the concentrated vorticity location. e above equation increases the vortex lift only when the leading-edge vortex (LEV) starts convection over the suction side.
By superposition of the loading, the total normal force coefficient is obtained as follows [14]: e chordwise force coefficient, pitching moment coefficient, and constants are considered according to Sheng guidelines [14,18].
where A is the rotor swept area, P describes the produced power from the product of the radius and the tangential force (rF t ), and F a is the axial force acting on the rotor.

Governing Equations.
e turbulent flow field around the wind turbine was computed by the finite volume method (FVM). e three-dimensional Reynolds-averaged Navier-Stokes (RANS) equations are formulated as follows [20]: e Reynolds stress term (− ρ � u i � u i ) was modeled based on the Boussinesq hypothesis according to equation (21), whereby the Reynolds stresses are concerned with the mean velocity gradients.
e realizable k-epsilon turbulence model (equations (22) and (23)) was used for its relatively low computational cost associated with the computation of the turbulent viscosity [21].
Regarding the temporal discretization, a second-order Crank-Nicolson method was applied for time integration. e PIMPLE algorithm was implemented for pressure-velocity coupling. is algorithm is a combination of SIMPLE and PISO algorithms. e matrices are solved by choosing the generalized geometric-algebraic multigrid (GAMG) for pressure and a preconditioned biconjugated gradient (PBiCG) for velocity [20]. e details of BEM and CFD data that are used in this paper for comparison are described in references [22,23].

Geometric Model.
In this study, the National Renewable Energy Laboratory (NREL) Phase VI reference wind turbine was modeled based on the sequence S settings represented in reference [24]. In this test sequence, an upwind, rigid turbine with a 0°cone angle was used, and the test was implemented for a wind speed, ranging from 5 m/s to 25 m/s. is case is comprised of a constant speed, constant pitch, and stallregulated wind turbine. It has a two-bladed rotor with twisted, tapered blades and shaped using the S809 airfoil [25,26]. e principal characteristics of the turbine are summarized in Table 1. e chord and twist distributions are illustrated in Figure 1. Figure 2 shows the cuboid computational domain, which consists of nonconformal hexahedral (cubic) grids, which allows local refinement in the areas of interest that significantly reduce the computational resources required for the simulation. Figure 2(a) expresses the boundary conditions for the wind turbine simulation. For the inlet condition, a uniform freestream velocity boundary was employed, whose magnitude was set to the examined velocity at the NREL experiment [24], and its direction was set parallel to the X-direction. e outlet condition was set as the pressure outlet boundary. A standard atmospheric pressure was considered for the reference pressure and imposed on the boundary. erefore, the dynamic pressure field is calculated considering this condition. No-slip boundary conditions were applied to the two lateral walls, the base, and the upper side of the domain.

Verification and Validation.
All the studied cases were simulated by considering a domain that is geometrically similar to the wind tunnel used in reference [24]. For this reason, the width and height of the domain section were set to 24.4 m and 36.6 m. e effect of upwind and downwind distances (d i , d o ) on the C P of the turbine was investigated as described in Table 2. is examination was performed at a wind speed equal to 7 m/s. While the mesh size kept the same, the domain length was extended according to Table 2. A comparison was made between the data and the longest domain (case 7). e results revealed that if the cases were longer than case 1, the difference was less than 0.1%. A compromise was reached between computational cost and accuracy, and thus case 2 (bolded) was chosen for this study. e effect of time step size was evaluated on the produced power coefficient (C P ). e results of the evaluation are given in Table 3. e required data are obtained for cases with 0°and 30°of yaw angles. It was performed to ensure involving phenomenon associated with yaw condition. For these tests, the wind speed was set to 7 m/s. e shortest time step with a value of Δθ � 0.216°was used as the reference case. In order to make data comparable, the ratio of C P /C P− ref is expressed in Figure 3, where the subscript "ref" refers to case 1. Figure 3 depicts for time steps smaller than Δθ � 0.432°; there was a difference of less than 2% in C P compared to the shortest time step size.
As shown in Figure 2, a structured hexahedral mesh was generated for the ALM-CFD calculation considering the mesh block structure. e computational domain was divided into four zones (see Figure 2(c)) to provide coarser grids in the far region. e cell size increases two times between each zone. e investigated cell size at the rotor disc is reported in Figure 4. As shown in this figure, a case with mesh Δx less than 0.125 m at the rotor disk can be chosen as the optimum mesh size with a balance between cost and accuracy. e independence of the following simulation results from the mentioned numerical factors was verified; the results were compared to the available data in the literature. In Figure 5, the numerical results for NREL Phase VI are reported together with the NREL experimental data [24]. e produced torque is reported as a function of wind speed, and as it can be observed, the computed aerodynamic torque agrees with the standard deviations for the tested wind speeds [27,28].     3.2. Yawed Case. Figure 6 shows the employed notation for yaw angle. Here, a positive yaw angle corresponds to a reduced inflow angle for the blade at 0°azimuth angle (12 o'clock). e right-handed Cartesian coordinate system is also depicted in this figure, whereby U ∞ is in the x-direction and tower is in the z-direction. In Figure 7(a), the produced aerodynamic torque of the turbine is compared for different yaw angles in U ∞ � 10 m/s. It is evident that an increase in the yaw results in the oscillation of produced power. It can be observed that by increasing the yaw angle to 10°, 20°, and 30°, the oscillation amplitude of produced torque increases 5%, 18%, and 20% of their average, respectively. is oscillation is addressed by the angle of attack variation that is due to the turbine yaw misalignment. Although the wind speed during rotation is almost constant, the yaw misalignment leads to oscillation in the angle between wind chord and wind flow oscillate during a rotor revolution.
is results in oscillation in both relative velocity and angle of attack. e magnitude of the geometrical relative velocity and the geometrical angle of attack are calculated using equations (24) and (25), respectively, as follows [29]: where c and φ denote the yaw angle and the azimuth angle and β stands for the pitch angle that includes the twist. e effect of turbine on the flow field is absent in the geometrical relative velocity and the geometrical angle of attack. Figure 7(b) illustrates how the oscillation of geometrical angle of attack varies with yaw angle during a rotor revolution. It should also be mentioned, while blade 1 is at its maximum angle of attack, the second blade is at its minimum.

Mathematical Problems in Engineering
An increment in yaw misalignment raises the oscillation amplitude while the average decreases.
is is a crucial phenomenon in wind turbine design. Nonetheless, simulation of it using a method like BEM is challenging at higher wind speed. e geometric angle of attack at the blade tip for the yaw angle of 30°is plotted in Figure 8(a) for different wind speeds. As shown in the figure, both the amplitude and the average of α geom grow with an increase in wind speed. e minimum angle of attack occurs at the azimuth angle of φ � 0°, after which it increases until its maximum at φ � 180°. e opposite is true when the yaw angle is in the negative direction. e geometrical relative velocity, V rel− geom , is also plotted in Figure 8(b). is figure shows that the variation of geometrical relative velocity varies inversely with the geometrical angle of attack, whereby the geometrical relative velocity reaches its maximum and minimum at azimuth angles of 0°and 180°, respectively.
As shown in Figure 8, turbines under yaw misalignment lead to a variation in the angle of attack. Hence, before employing the Sheng DS model in the study of yawed cases, an assessment was established between the Sheng DS model results and the experimental measurements reported in [14]. For this part, the simulation was performed for an isolated airfoil, whereas in the rest of this paper, an entire turbine has been simulated using the ALM. e results of the ramp-up test at a reduced pitch rate (r � _ αc/2U ∞ ) of 0.005 are plotted in Figure 9(a). e results for a single S809 produced by the employed DS model were very close to the experimental data in terms of the normal force. It predicted the stall onset just 2°earlier, and the effect of the vortex emitted at 27°was less than that in the experimental data; however, the results were still in good agreement. Furthermore, the results of a single oscillatory airfoil were compared with the experimental data, as reported in Figure 9(b). e hysteresis loop is obtained for an oscillatory airfoil with α � 15 ± 10°and a reduced frequency of κ � ωc/2U ∞ � 0.074. Again, the predicted stall was slightly earlier. It should be noted that the employed DS model was initially developed for helicopter applications, and even the modified version by Sheng inherited the main characteristics of the original model, e.g., employing the Kirchhoff equation, which is more suitable for thinner airfoils. Despite this, employing the Sheng DS model established significant capability in modeling delays and downstroke movement in the dynamic stall.
In Figure 10, the aerodynamic torque calculated using the current ALM code is compared with UBEM and CFD results, the details of which are available in [22,23]. Data are presented for different wind speeds of 10 m/s, 13 m/s, and 15 m/s corresponding to TSR equal to 3.8, 2.9, and 2.5, respectively. All three figures were obtained for the NREL Phase VI wind turbine under 30°of yaw angle. In general, the accuracy and reliability of the presented ALM method were confirmed by CFD computation as the percentage error in the average aerodynamic torque between CFD, and the ALM calculation tools were very limited. e results can be compared in terms of amplitude and phase of the oscillating torque for further elaboration.
As can be observed, UBEM underestimated the torque for U ∞ � 10 m/s and overestimated it for U ∞ � 15 m/s; nonetheless, there was an outstanding improvement in the amplitude by ALM. It should be mentioned that UBEM uses yaw models to consider the nonorthogonality of flow to the rotor plane [30], whereas the ALM model obtains the inflow velocity from the resolved flow field; thus, all flow angles are inherently concerned.
Regarding the phase of oscillating torque, a phase shift error of about 30°and 15°appeared between UBEM and CFD results at 10 m/s and 13 m/s, respectively. Nevertheless, as the wind speed increased, the phase error decreased towards zero. Considering the angle of attack, the blades are subjected to dynamic stall at a wind speed of 10 m/s in some portions of a revolution. As the wind speed increases to 13 m/s, the blade is subjected to dynamic stall almost during the entire revolution. An outstanding improvement in decreasing the phase shift error appeared in the ALM results compared to UBEM. e phase was predicted accurately for U ∞ � 10 m/s with a partial dynamic stall. However, the accuracy decreased slightly for U ∞ � 13 m/s as there was an incremental oscillation after φ � 90°, caused by a delay in predicting the dynamic stall of the second blade. As the wind speed increased to 15 m/s, the high angle of attack during the revolution led to a complete flow separation and the formation of TEV. While the ALM uses the tabulated airfoil data, this three-dimensional phenomenon affects the inflow, resulting in a slightly inaccurate prediction by ALM. e results obtained for the NREL Phase VI reference wind turbine with different static yaw misalignments are presented in Figure 11. Cases with yaw angles ranging from 0°to 30°were simulated in both negative and positive yaw directions. Square symbols show the average power, and the minimum and maximum values over the revolution are indicated by error bars. As shown, the results for negative and positive yaw angles are relatively symmetrical. e results indicate an average power reduction of 25% due to 30°yaw misalignment. e variation of power with yaw angle P(c) was suggested to be approximated by P(0°)cos c χ [31]. e exponent χ was suggested to be equal to 3; however, it is only valid if the axial induction distribution is small compared to U ∞ cos c [31]. It was stated that measurements including wind tunnel models and full-scale turbines have indicated the exponent χ could vary between 1.88 and 5.14 [32]. In the current case, χ � 2 shows to be in good agreement for the studied yaw misalignments by the ALM results, although it is not entirely fit.
Taking the amplitude of oscillations into consideration, it shows to grow continuously until 30°, after which there is no clear trend. is is related to the angles of attack as for cases with yaw angles greater than 30°, the blade is under stall condition in a portion of its revolution. e blades being under stall conditions lead to instabilities for further investigations. In other words, the amplitudes for these conditions are case-dependent, as they is a result of two stalled blades. Figure 12 illustrates the wake structure as viewed from the side (x-z plane). e isosurface of Q visualizes the instantaneous vorticity, where the Q criterion is defined as follows [33]:

Mathematical Problems in Engineering
where Ω and S represent the symmetric and antisymmetric parts of the velocity gradient tensor ∇u, respectively: e flow field around a flow-aligned turbine and a 30°y awed case at a wind speed of 10 m/s is shown in Figures 12  and 13, respectively. For the flow-aligned turbine, the wake was initially developed behind the rotor blades in the form of a well-defined helical geometry. e wake structure was conserved until four rotor diameters downstream and then evolved into the turbulent wake. e strength of the helical wake vorticities was azimuthally symmetric; hence, the contour in Figure 12(b) was illustrated only for one time step.
For the yawed case, the contour at the blade section is illustrated for each 30°of revolution as indicated in Figure 13(b) I-XII . e comparison of yawed cases with the flow-aligned case indicates that the yaw misalignment yielded significantly different wake evolution. A periodic variation in the wake strength with respect to the azimuthal angle was manifested. e periodical deformation of the wake vorticities strength is addressed by the variation of AOA and relative velocity during revolutions. e corresponding AOA and relative velocity can be observed in Figure 8. As shown in Figure 13(b) I , a more robust wake was emitted when the blades were at the azimuth position of 0°. is position corresponds to the highest relative velocity during a revolution.
From a downstream wind turbine perspective, axial velocity distribution directly influences the inflow condition of the downstream wind turbines. e inflow condition has a crucial role in determining the aerodynamic and fatigue loads of the downstream wind turbine [34]. In Figures 14  and 15 , the streamwise velocity contours behind the aligned turbine are compared with the yawed wind turbine to ascertain the influence of emitted wake structure on the downstream velocity field. e data were reported at x/D � 1, 4, and 8 as indicated in Figures 14(a) and 15(a), where x is the downstream distance and D is the rotor diameter. Simulation indicated that the propagation of the wake lasted longer in the case of the aligned rotor. As it can be observed in contours at x/D � 8, the portion of contour in dark blue is broader for the aligned case. e wake evolved in the form of a helical geometry with similar strength for the aligned case, as shown in Figure 14(a).
us, the distribution of velocity deficit is axisymmetric concerning the wake age. As the wind turbine underwent 30°of yaw angle, a highly skewed deformed wake structure was developed (see Figure 15). According to Figure 15(b), the axial velocity distribution was no longer axisymmetric concerning the wake age because of the   periodically oscillating wake vortices due to oscillating AOA and relative velocity. As strong wake vortices pass the downstream region, a high induced velocity region appeared at the azimuth angle of 240°. is agrees with the observation of Bastankhah and Porté-Agel [35] where the wake centre, defined as the point where the velocity deficit is maximum at each downwind location, was changed according to yaw angle. us, the downstream wind turbine receives significantly unsteady and asymmetric inflow conditions and blades suffer from a time-dependent aerodynamic load.

Conclusions
An actuator line model was implemented in OpenFOAM flow solver. e Shen's tip treatment to consider the tip effect and Sheng's LB dynamic stall model to consider the dynamic stall effects were applied. e ALM simulation results for the aligned case were in excellent agreement with available experiment, and the overall trend was reproduced for the output torque with respect to wind speed. An increment in yaw angle leads to blade advancing and retreating effect whereby the angle of attack oscillates. e variation of produced power was investigated in the following three conditions: (1) e turbine at a wind speed of 10 m/s and under 30°o f yaw represents an understall condition. In this condition, ALM results were in general matching in terms of both phase and amplitude.
(2) e turbine at a wind speed of 13 m/s and under 30°o f yaw is a condition that the blade is subjected to dynamic stall almost during the entire revolution. Employing the Sheng DS model established a significant capability in modeling delays and downstroke movement in the dynamic stall.
(3) e turbine at a wind speed of 15 m/s and under 30°o f yaw represents a deep stall-dominated flow. e phase was predicted accurately with a minor overestimation. Considering the fact that ALM reads airfoil coefficients from tabulated data, further modification in the model is required if this stalldominated condition is in the scope of a study.
In general, despite the small number of cells, the ALM code is capable of predicting produced power with satisfying accuracy before deep stall-dominated conditions. Employing the DS model improved the prediction of phase delay in power variation.
Providing knowledge about turbine wakes is an advantage of ALM over BEM, which is a priority in wind farm design. According to the results, even using 400k cells, the ALM code could fairly simulate the skewed helical geometry of wake and capture the effect of advancing and retreating blade on the wake. It was found that the dominated flow by the helical wake lasts longer in the aligned case. On the other side, the deficit velocity is more concentrated for the yawed case. Considering these points can be exploited for wind farm optimization, ALM potentially promises to be utilized with less computational cost than CFD. Normal force coefficient for the trailing-edge separation (-) C ] N : Normal force coefficient for the vortex shedding (-) C P :

Nomenclature:
Power coefficient (-) C 1 : Realizable k-epsilon turbulence model coefficient (-) C 2 : Realizable k-epsilon turbulence model coefficient (-) D: Effective diffusivity (-) Deficiency function for the impulsive force coefficient (-) g: Gravity acceleration (m/s 2 ) G: Turbulent kinetic energy production rate due to the anisotropic part of the Reynolds stress tensor (m 2 ·s − 3 ) M: Mach number (-) p: Pressure (Pa) r: Distance from actuator point to the point where the force is applied (m) Section width (m).
Greek: α: Local angle of attack _ α: Time derivative of angle of attack (rad/s) α ′ : Delayed angle of attack (rad) α E : Equivalent angle of attack (rad) α cr : Critical stall-onset angle of attack (rad) α ds 0 : Constant critical stall-onset angle of attack (rad) α ss : Static stall-onset angle (rad) α 1 : Angle of attack for the breakpoint of the flow separation (rad) β: Pitch angle c: Yaw angle Δ: Step change in angle of attack or in time ε: Projection width (m) ϵ: Turbulent kinetic energy dissipation rate (m 2 ·s − 3 )