An extended k − ε model for wake-flow simulation of wind farms

The Reynolds-averaged Navier-Stokes approach coupled with the standard k − ε model is widely utilized for wind-energy applications. However, it has been shown that the standard k − ε model overestimates the turbulence intensity in the wake region and, consequently, overpredicts the power output of the waked turbines. This study focuses on the development of an extended k − ε model by incorporating an additional term in the turbulent kinetic energy equation. This term accounts for the influence of turbine-induced forces, and its formulation is derived through an analytical approach. To assess the e ff ectiveness of the proposed model, we begin by analyzing the evolution of normalized velocity deficit and turbulence intensity in the wake region, and the normalized power of the waked turbines. This investigation involves a comparison of the predictions against results from large-eddy simulations in three validation cases with di ff erent layouts. We then simulate a wind farm consisting of 30 wind turbines and conduct a comparative analysis between the model-predicted normalized streamwise velocity and wind-tunnel measurements. Finally, to conclude our assessment of the proposed model, we apply it to the operational wind farm of Horns Rev 1 and evaluate the obtained normalized power with the results from large-eddy simulations. The comparisons and validations conducted in this study prove the superior performance of the extended k − ε model compared to the standard version.


Introduction
Harvesting renewable energy from the atmospheric boundary layer (ABL) through wind turbines is one of the most promising methods for clean power production.Following the policies for the mitigation of global warming, wind will be the prominent power source by 2050, generating one-third of the total electricity demand.This translates into scaling up the annual wind-energy installations by four times in the upcoming decade [1][2][3].Such a deployment requires understanding and accurate prediction of the multi-scale fully-coupled interactions between the ABL and wind farms.Within wind farms, the turbulent flows -known as the wakes -that form behind turbines are associated with a strong velocity deficit and elevated turbulence level compared to the upstream flow.The velocity deficit decreases with distance owing to the turbulent mixing mechanism.However, due to the short inter-turbine spacing in wind farms, the velocity does not completely recover to the upstream value.This results in a power loss of the waked turbines and an increase in their structural loads (see the review of Refs.[4][5][6] and references therein).Therefore, accurate predictions of turbine wakes and power losses are essential for the design of new wind farms and the optimization of the annual energy production and loads.
Owing to the growth of computing power, computational fluid dynamics (CFD) is broadly used for wind-farm modeling.CFD techniques, i.e., Reynolds-averaged Navier-Stokes (RANS) models [7] and large-eddy simulation (LES) [8], can capture the physics of flow within wind farms to a much higher extent compared to the simple analytical-empirical wake models [9][10][11].LES provides a high level of physical fidelity and compares well with data from experiments and field measurements [12,13], nevertheless, it is costly and more suitable for problems with where u i , p, S i j , R i j , and f i denote the mean velocity, mean pressure, mean rate-of-strain, the Reynolds stress, and turbine-induced forces, respectively.The fluid's density and kinematic viscosity are represented with ρ and ν, respectively.The Reynolds-stress term is formulated through the eddy-viscosity hypothesis, i.e., R i j = −2ν T S i j + 2/3kδ i j , where, ν T and δ i j are the eddy viscosity and the Kronecker delta, respectively [23,24].
To model ν T , an ensemble of eddy-viscosity models is available in the literature (see, e.g., Refs.[25][26][27][28][29]).Here, we focus on the widely used standard k − ε model in which the eddy viscosity is formulated as ν T = C µ k 2 /ε, where C µ is a model constant.The transport equations of k in our extended k − ε model is where σ k is a constant.The additional term (S k ) which is absent in the standard k − ε model, accounts for the impact of turbine forces that needs to be parameterized.

Parametrization of the turbine-induced terms in the RANS framework
We use the standard actuator-disk model without rotation to calculate the turbine-induced forces ( f i ) based on the disk-averaged velocity (u d,i ) by using the disk-based thrust coefficient C ′ T = C T /(1 − a) 2 , in which a is the turbine's induction factor [17,30].The formula to calculate f i is where n i denotes the unit vector perpendicular to the disk.The distribution of the forces is proportional to the frontal surface of the computational mesh (A cell ) that falls within the rotor area.Here, γ j,l is the fraction of the area overlap between the cell at a grid point ( j, l) and the rotor.It has a value of 1 and 0 for all cells with and without overlap with the rotor, respectively.For cells having a partial overlap with the rotor, γ j,l is equal to the fraction of area overlap [30].
The transport equation of k is derived through a product of fluctuating velocity by the subtraction of the Reynoldsaveraged momentum equation from the instantaneous momentum equation, and then taking the average of the resulting equation.For a yaw-aligned rotor, assuming that the free stream wind is aligned with the x-direction, a sink term corresponding to the turbine-induced forces appears in the k equation as By assuming u ′2 d,x ≈ 2/3k d and u ′3 d,x ≈ u ′2 d,x 3/2 [23], S k can be further simplified as where k d is the disk-averaged TKE.

Description of validation cases and numerical setups
Here, we define three validation cases to test the performance of the extended k − ε model by comparing its predictions against those from our LESs and the standard k − ε model.The validation cases have six turbines with a rotor diameter, hub height, and induction factor of 80 m, 70 m, and 0.25, respectively.The inflow velocity and inflow turbulence intensity at the hub height are 8 m/s and 5.8%, respectively.In case 1, shown in Figure 1(a), the inter-turbine spacing is 7d o which is a typical value for wind-farm layouts [31].To include a case with more severe wake effects compared to case 1, we define case 2 with an inter-turbine spacing of 5d o (Figure 1(b)).To include a case with partial wakes, the turbines on the even rows are relocated in the y-direction by d 0 in case 3 (Figure 1(c)).
As our RANS framework, we use the simpleFoam solver, available in OpenFOAM-v2112 [32].We implement the new turbulence model, to solve the above-mentioned governing equations and simulate the turbine wakes and power losses under neutral ABL.The computational domains have a size of 4400 m × 400 m × 355 m in the x, y, and z directions, respectively.In all three cases, the first turbine is located at a distance of 5d 0 from the domain's inlet.
Our grid-convergence study shows a resolution of 234 × 40 × 58 in the x, y, and z direction is sufficient for all the simulations.We enhance the mesh resolution in the x direction at the turbines' location.We also maintain a uniform division of the domain in the spanwise direction while employing a stretching grid in the z direction.We apply the cyclic boundary conditions at the side boundaries.A pressure-outlet condition for the pressure, and a zero-gradient boundary condition for all other quantities are utilized at the outlet.At the inlet and top boundaries, for the velocity and turbulence quantities, values corresponding to a logarithmic boundary-layer flow [31] are utilized, and for pressure, a zero-gradient condition is used.A rough wall boundary condition with wall functions for the turbulence quantities is imposed at the bottom wall [32].
The LES framework and the numerical setups utilized for the simulation of the validation cases are similar to our previous works, and for the sake of brevity, we do not include further details on those.Readers are referred to Refs.[32,33], for more information.The LES framework has been validated and extensively used in prior studies, e.g., see Refs. [34][35][36][37][38].

Results and discussion
In this section, first, a thorough analysis is performed on the performance of the extended k − ε model when applied to the validation cases.The predicted normalized velocity deficit (∆u/u h = (u in − u)/u h with u in and u h as the inflow velocity and the hub-height inflow velocity, respectively) and turbulence intensity (I = √ 2k/3/u h ) are compared against those from LES.Then, the vertical profiles of normalized streamwise velocity (u/u h ) downstream of several rows in a 10×3 array of turbines are compared with the wind-tunnel measurements performed by Chamorro and Porté-Agel [39].Finally, we apply our model to the HR1 wind farm and evaluate its predictions against LES data from Wu and Porté-agel [40].The extended k − ε model would be considered successful if it could outperform the standard k − ε in all the mentioned cases.power losses under neutral ABL.The computational domains have a size of 4400 m × 400 m × 355 m in the x, y, and z directions, respectively.In all three cases, the first turbine is located at a distance of 5d 0 from the domain's inlet.
Our grid-convergence study shows a resolution of 234 × 40 × 58 in the x, y, and z direction is sufficient for all the simulations.We enhance the mesh resolution in the x direction at the turbines' location.We also maintain a uniform division of the domain in the spanwise direction while employing a stretching grid in the z direction.We apply the cyclic boundary conditions at the side boundaries.A pressure-outlet condition for the pressure, and a zero-gradient boundary condition for all other quantities are utilized at the outlet.At the inlet and top boundaries, for the velocity and turbulence quantities, values corresponding to a logarithmic boundary-layer flow [31] are utilized, and for pressure, a zero-gradient condition is used.A rough wall boundary condition with wall functions for the turbulence quantities is imposed at the bottom wall [32].
The LES framework and the numerical setups utilized for the simulation of the validation cases are similar to our previous works, and for the sake of brevity, we do not include further details on those.Readers are referred to Refs.

Results and discussion
In this section, first, a thorough analysis is performed on the performance of the extended k − ε model when applied to the validation cases.The predicted normalized velocity deficit (∆u/u h = (u in − u)/u h with u in and u h as the inflow velocity and the hub-height inflow velocity, respectively) and turbulence intensity (I = √ 2k/3/u h ) are compared against those from LES.Then, the vertical profiles of normalized streamwise velocity (u/u h ) downstream of several rows in a 10×3 array of turbines are compared with the wind-tunnel measurements performed by Chamorro and Porté-Agel [39].Finally, we apply our model to the HR1 wind farm and evaluate its predictions against LES data from Wu and Porté-agel [40].The extended k − ε model would be considered successful if it could outperform the standard k − ε in all the mentioned cases.power losses under neutral ABL.The computational domains have a size of 4400 m × 400 m × 355 m in the x, y, and z directions, respectively.In all three cases, the first turbine is located at a distance of 5d 0 from the domain's inlet.
Our grid-convergence study shows a resolution of 234 × 40 × 58 in the x, y, and z directions is sufficient for all the simulations.We enhance the mesh resolution in the x-direction at the turbines' location.We also maintain a uniform division of the domain in the spanwise direction while employing a stretching grid in the z-direction.We apply the cyclic boundary conditions at the side boundaries.A pressure-outlet condition for the pressure, and a zero-gradient boundary condition for all other quantities are utilized at the outlet.At the inlet and top boundaries, for the velocity and turbulence quantities, values corresponding to a logarithmic boundary-layer flow [33] are utilized, and for pressure, a zero-gradient condition is used.A rough wall boundary condition with wall functions for the turbulence quantities is imposed at the bottom wall [34].
The LES framework and the numerical setups utilized for the simulation of the validation cases are similar to our previous works, and for the sake of brevity, we do not include further details on those.Readers are referred to Refs.[34,35], for more information.The LES framework has been validated and extensively used in prior studies, e.g., see Refs.[36][37][38][39][40].

Results and discussion
In this section, first, a thorough analysis is conducted on the performance of the extended k − ε model when applied to the validation cases.The predicted normalized velocity deficit (∆u/u h = (u in − u)/u h with u in and u h as the inflow velocity and the hub-height inflow velocity, respectively) and turbulence intensity (I = √ 2k/3/u h ) are compared against those from LES. Aiming to provide a more comprehensive and quantifiable understanding of the physical phenomena under consideration, we also examine the evolution of the ratio of eddy viscosity to the inflow value (ν T /ν T 0 ) in case 1, predicted through the standard and extended k − ε models.Afterward, the vertical profiles of normalized streamwise velocity (u/u h ) downstream of several rows in a 10 × 3 array of turbines are compared with the wind-tunnel measurements performed by Chamorro and Porté-Agel [41].Finally, we apply our model to the HR1 wind farm and evaluate its predictions against LES data from Wu and Porté-agel [42].The extended k − ε model would be considered successful if it could outperform the standard k − ε in all the mentioned cases.

Evaluating the extended k − ε model on the validation cases
In Figure 2, two-dimensional representations are adopted to show the normalized velocity deficit at the hub height in case 1, obtained from LES, the standard k−ε model, and the extended k−ε model, from top to bottom.Compared to the results from LES and the extended k − ε model, the standard k − ε model underestimates the velocity deficit behind all turbines, with the largest deviation observed for the most upstream one.It is essential to recognize that turbines exert a substantial pressure gradient upon the flow within their very vicinity.This circumstance poses a significant challenge to the eddy-viscosity models due to their inherent limitation in capturing pressure-velocity correlations [14].Consequently, this limitation impacts the accuracy of predictions rendered by both the standard k − ε model and the extended k −ε model in the near wake.In this context, focusing on the downstream distribution of normalized velocity deficit, i.e., at x/d 0 ≥ 5, shows that the extended k − ε model has successfully captured a distribution close to LES which is essential for an accurate prediction of normalized power of the waked turbines.
To shed light on the phenomena observed in the previous figure, and to comprehend the physics of flow in turbine wakes as captured by the standard and extended k − ε models, we study the evolution of eddy viscosity in case 1. Figures 3(a 1 ) and (a 2 ) illustrate the hub-height level contours of the ratio of eddy viscosity to the inflow value (ν T /ν T 0 ).As seen in these figures, the term introduced in the TKE equation successfully addresses the well-known issue associated with the standard k − ε model regarding the overestimation of the eddy viscosity and mixing [7,17].Additionally, Figure 3(b) presents the rotor-averaged ν T /ν T 0 , where rotor-averaging is performed over the intersection area of the x-normal planes and a hypothetical cylinder starting from the inlet with the same diameter as that of the turbines with its axis at the hub height.It is evident that the extended k − ε model significantly reduces ν T /ν T 0 in the wake region at all x-positions.These findings are reasonable to extrapolate to cases 2 and 3 and, for the sake of brevity, we omit the results for these cases.) fails to capture the double peak structure of the turbulence intensity behind the turbines.This is consistent with the distribution of eddy viscosity presented in Figure 3(a 1 ), where the standard k − ε model exhibited a single-peak behavior with a maximum at the center of the wake across all locations.In contrast, the extended k − ε model (Figures 3(a 2 )) demonstrated an opposite behavior with a double-peak structure, leading to a reduction of ν T /ν T 0 in the central part of the wake compared to the sides.Owing to the overestimation of the eddy viscosity and enhanced mixing, the standard k − ε model predicts a strong development of turbulence intensity in turbine wakes, especially behind the first turbine.Turning to the extended k − ε model, the new term included in the TKE equation acts as a sink at the location of the turbines and plays an important role in controlling the production of turbulent kinetic energy in the wake region.
To assess the performance of the extended k − ε model more quantitatively, the normalized power of turbines is analyzed.The success level of a model in the prediction of the normalized power is highly dependent on the accuracy of the obtained disk-averaged velocity.Therefore, with the results of the previous figures in mind, we expect the extended k − ε model to outperform the standard version.Figure 5 compares the turbines' normalized power predicted by the extended k − ε model with the values obtained from LES and the standard k − ε model.Apart from the second turbine, the extended k − ε model shows a good agreement with the LES and outperforms the standard k − ε model for all the waked turbines.The overprediction of the normalized power for the second turbine can be primarily attributed to a substantial velocity recovery within the wake of the most upstream turbine, which results from the strong turbulent mixing predicted through RANS models.
Figures 6 and 7 depict the contour plots of normalized velocity deficit and turbulence intensity at the hub height in case 2. As mentioned earlier, this case has a stronger wake effect compared to case 1 due to a shorter inter-turbine spacing.The predicted normalized velocity deficit by the extended k − ε model (Figure 6(c)) is closer to the LES results (Figure 6(a)), especially for the waked turbines.This shows that the improvement compared to the standard k − ε model would be higher in cases with strong wake effects.Similarly, the term introduced in the TKE equation is helping to hamper the overprediction of the turbulence intensity behind the turbines.Also, the double-peak structure of the turbulence intensity behind the turbines is captured by the extended k − ε model (Figure 7(c)) similar to LES (Figure 7(a)).
The normalized power of turbines in case 2, obtained from LES, the extended k − ε model, and its standard version are shown in Figure 8. Due to the higher wake losses of case 2 compared to case 1, the normalized powers of all                  waked turbines are degraded in this case.A good agreement between the results of LES and the extended k − ε can be observed, especially for turbines 3 to 6.Much like the scenario witnessed in case 1, the strong development of RANS-predicted turbulence intensity behind the very first turbine due to overestimation of eddy viscosity and, as a result, a fast recovery of velocity, underlies the higher error of the normalized power for the second turbine.
Figures 9 and 10 illustrate the hub-height contour plots of normalized velocity deficit and turbulence intensity in case 3.In this case, downstream-located turbines experience partial wakes and mild full wakes.Thus, as observed in Figure 9, the normalized velocity deficit is smaller than that in the cases investigated earlier.In terms of matching with LES results, in a similar manner as the previous cases, we can observe a good agreement for the extended k − ε model.According to Figure 10, and by focusing on the LES results, one can see that the added turbulence in the wake region of case 3 is also smaller compared to cases 1 and 2. The overestimation of the turbulence intensity by the standard k − ε model can be observed in this figure.On the contrary, the turbulence intensity obtained from the extended k − ε model shows a good agreement with the LES results.
Figure 11 illustrates the normalized power in case 3 calculated through LES, the standard k − ε model, and our extended version.Owing to the fact that turbine 2 is only affected by a partial wake from turbine 1, its normalized      The normalized power of turbines in case 2, obtained from LES, the extended k − ε model, and its standard version are shown in Figure 7. Due to the higher wake losses of case 2 compared to case 1, the normalized powers of all waked turbines are degraded in this case.A good agreement between the results of LES and the extended k − ε can be observed, especially for turbines 3 to 6.Much like the scenario witnessed in case 1, the stronger development of RANS-predicted turbulence intensity behind the very first turbine underlies the higher error of the normalized power for the second turbine.
Figures 8 and 9 illustrate the hub-height contour plots of normalized velocity deficit and turbulence intensity in 6 (Figure 7(a)).power is significantly higher than the second turbines in cases 1 and 2, and a perfect match between the three models is observed.Focusing on the other waked turbines, the extended k − ε model performs satisfactorily, with a good agreement with the results from LES.
Table 1 provides the average value of the normalized power of waked turbines (No. 2 -6) calculated through the three models of LES, standard, and extended k − ε.The maximum errors compared to LES occur in case 2 with the most severe wake effect for both RANS models, however, the extended k − ε model has a considerably lower error compared to the standard version.The reduction in error is maximum (37.50%) for case 2, and in case 3 where the           case 3.In this case, downstream-located turbines experience partial wakes and mild full wakes.Thus, as observed in Figure 8, the normalized velocity deficit is smaller than that in the cases investigated earlier.In terms of matching with LES results, in a similar manner as the previous cases, we can observe a good agreement for the extended k − ε model.
According to Figure 9, and by focusing on the LES results, one can see that the added turbulence in the wake region of case 3 is also smaller compared to cases 1 and 2. The overestimation of the turbulence intensity by the standard k − ε model can be observed in this figure.On the contrary, the turbulence intensity obtained from the extended k − ε model shows a good agreement with the LES results.power is significantly higher than the second turbines in cases 1 and 2, and a perfect match between the three models 178 standard k − ε model performs well, the reduction in error is 8.95%.Therefore, we can infer that the strategy for including the analytically-derived term in the TKE equation is successful and the extended model performs well when applied to an ensemble of wind-farm cases with different layouts and levels of wake strength.

Evaluating the extended k − ε model against wind-tunnel experiments
Building upon the promising results demonstrated by the extended k − ε model when applied to the validation cases, we now provide further exploration of the model's capabilities.To this end, the wind-tunnel experiment under neutral stratification performed by Chamorro and Porté-Agel [41] is selected as the reference.Their experimental    Table 1 provides the average value of the normalized power of waked turbines (No. 2 -6) calculated through the 177 three models of LES, standard, and extended k − ε.The maximum errors compared to LES occur in case 2 with the 178 most severe wake effect for both RANS models, however, the extended k − ε model has a considerably lower error 179 compared to the standard version.The reduction in error is maximum (37.50%) for case 2, and in case 3 where the 180 standard k − ε model performs well, the reduction in error is 8.95%.Therefore, we can infer that the strategy for 181 including the analytically-derived term in the TKE equation is successful and the extended model performs well when 182 applied to an ensemble of wind-farm cases with different layouts and levels of wake strength.Table 1 provides the average value of the normalized power of waked turbines (No. 2 -6) calculated through the three models of LES, standard, and extended k − ε.The maximum errors compared to LES occur in case 2 with the most severe wake effect for both RANS models, however, the extended k − ε model has a considerably lower error compared to the standard version.The reduction in error is maximum (37.50%) for case 2, and in case 3 where the standard k − ε model performs well, the reduction in error is 8.95%.Therefore, we can infer that the strategy for including the analytically-derived term in the TKE equation is successful and the extended model performs well when applied to an ensemble of wind-farm cases with different layouts and levels of wake strength.setup consisted of a 10 × 3 array of miniature wind turbines with a hub height of 0.83d 0 , a streamwise spacing of 5d 0 , and spanwise spacing of 4d 0 , where d 0 = 0.15 m.The reported roughness height and boundary-layer depth in the experiments were 19.98 × 10 −5 d 0 and 4.5d 0 , respectively [41,43].To simulate this wind farm with our RANS framework, a computational domain with a size of 95d 0 × 12d 0 × 4.5d 0 in the streamwise, spanwise, and vertical directions, respectively, is created.The layout of the wind farm is schematically shown in Figure 12(a).Similar to the validation cases, the first turbine row is placed at a distance of 5d 0 from the domain inlet.A grid with a resolution of 450 × 96 × 58 is used in the streamwise, spanwise, and vertical directions, respectively.We refine the mesh at the turbines' location in the x-direction and the domain is uniformly divided in the y-direction.Additionally, we utilize a stretching grid in the z-direction.Our assessment of grid independence demonstrates that this resolution effectively captures the key features of turbine wakes.The inflow characteristics match those from the experiment, and the boundary conditions are similar to the ones described in section 2.2.As depicted in Figure 12(b), to simulate the wind farm, we use the values of C ′ T for the different turbine rows reported in Ref. [43]. Figure 13 shows the vertical distributions of the normalized streamwise velocity obtained from the standard and extended k − ε models compared against wind-tunnel measurement data [41].The profiles are extracted at a distance validation cases, the first turbine row is placed at a distance of 5d 0 from the domain inlet.A grid with a resolution of 450 × 96 × 58 is used in the streamwise, spanwise, and vertical directions, respectively.We refine the mesh at the turbines' location in the x direction and the domain is uniformly divided in the y direction.Additionally, we utilize a stretching grid in the z direction.Our assessment of grid independence demonstrates that this resolution effectively captures the key features of turbine wakes.The inflow characteristics match those from the experiment, and the boundary conditions are similar to the ones described in section 2.     of 3d 0 behind the turbine rows 1-7 and 10 [43].Focusing on the predictions of the standard k − ε model, one can see that it generally overestimates the normalized streamwise velocity in the area from bottom-to top-tip heights, and this agrees with the earlier observations.In contrast to the standard k − ε model, there is a satisfactory agreement between the results of the extended k − ε and measurement data within the bottom-tip to the top-tip range.This observation is also in agreement with previous findings showing that the extended k − ε model predicts more accurate disk-averaged velocities and normalized power of turbines.

Applying the extended k − ε model to the Horns Rev 1 (HR1) wind farm
To evaluate the performance of the extended k − ε model on operational wind farms, we apply it to the HR1 wind farm and compare the obtained normalized power with LES data from Ref. [40] for a hub-height inflow velocity of 8 m/s and ambient turbulence intensity of 7.7%.The HR1 wind farm holds eighty Vestas V-80 2 MW turbines with a hub height of 70 m and a rotor diameter of 80 m.The spacing between the consecutive turbines for wind direction of 270 • , 222 • , and 312 • are 7d 0 , 9.3d 0 , and 10.4d 0 , respectively.Therefore, among these full-wake conditions, the most severe wake effect would be observed for a wind direction of 270 • .We simulate the entire wind farm with our RANS framework coupled with the standard and extended k − ε models.For the sake of brevity, we omit more information on the computational domain and boundary conditions.Readers may refer to Ref. [42], for more information on the methodology adopted for the simulations.

Conclusions
The present study aimed to propose an extended k − ε model for the simulation of wakes and power losses in wind farms.It is known that the inherent shortcomings associated with the widely-utilized standard k − ε model lead to overestimation of eddy viscosity and turbulence intensity in the wake region and, as a result, an overprediction of the normalized power of the waked turbines.To solve this issue, we included an additional term resembling the impact of the turbine forces in the transport equation of turbulent kinetic energy.The new term was derived based on physical principles through an analytical approach.To assess the performance of the proposed model, initially, we studied three validation cases with different layouts under full and partial wake conditions and compared the quantities of interest, i.e., normalized velocity deficit, turbulence intensity, and turbines' normalized power, predicted by the standard and extended k −ε models against LESs performed in-house [34,35].Then, we modeled a wind farm consisting of a 10×3 array of turbines and evaluated the success of the extended k − ε model in predicting the streamwise velocity in turbine wakes through comparison with wind-tunnel measurement from Ref. [41].Finally, we examined the performance of our model by applying it to the operational wind farm of HR1 and compared its predictions with LES data from Ref. [42].
The investigations showed that the incorporation of the term associated with the turbine-induced forces in the transport equation of the turbulent kinetic energy managed to indirectly reduce the eddy viscosity and, consequently, the mixing, in turbine wakes, as a well-known challenge linked to the standard k − ε model.The turbulence level behind the turbines, predicted by the extended k − ε, compared well with the LES results, and unlike the standard k − ε model, our extended model could capture the double peak structure of the turbulence intensity behind the turbines.As a result, the extended model could solve the issue with the standard k − ε model in terms of the underestimation of velocity deficit in turbine wakes.Then, we focused on the normalized powers obtained from the standard and extended versions of the k − ε model and compared them with LES.In all three validation cases, the extended k − ε model outperformed the standard k − ε.In the context of comparison against wind-tunnel measurements, it became evident that the normalized streamwise velocity in the wake region obtained from the extended model exhibited satisfactory agreement with measurement data across the vertical distribution from bottom-to top-tip heights.When applied to the operational wind farm of HR1, the relative error between the predictions of the extended k − ε and the LES data [42] were below 8% for three wind directions leading to the full-wake condition.The comparisons and validations conducted in this study proved the high potential of the proposed model for utilization by the wind-energy community.Applications of the present model to other complex turbulent flow problems (e.g.hilly terrain and thermally stratified ABL) could be considered as future work to further enrich the assessment of the proposed model.

Figure 1 :
Figure 1: Schematics of the validation wind-farm cases: (a) Case 1 with a 7d 0 streamwise spacing, (b) case 2 with a 5d 0 streamwise spacing, and (c) case 3 with a 7d 0 streamwise spacing where even turbines are relocated by d 0 in the y-direction.Here, the black rectangles show the turbines' location.Inlet and outlet boundaries of the domain are located at x/d 0 = −5 and 50, respectively.The black dashed lines demonstrate the cyclic boundaries.

3. 1 . 4 Figure 1 :
Figure 1: Schematics of the validation wind-farm cases: (a) Case 1 with a 7d 0 streamwise spacing, (b) case 2 with a 5d 0 streamwise spacing, and (c) case 3 with a 7d 0 streamwise spacing where even turbines are relocated by d 0 in the y-direction.Here, the black rectangles show the turbines' location.Inlet and outlet boundaries of the domain are located at x/d 0 = −5 and 50, respectively.The black dashed lines demonstrate the cyclic boundaries.

3. 1 . 4 Figure 1 :
Figure 1: Schematics of the validation wind-farm cases: (a) Case 1 with a 7d 0 streamwise spacing, (b) case 2 with a 5d 0 streamwise spacing, and (c) case 3 with a 7d 0 streamwise spacing where even turbines are relocated by d 0 in the y-direction.Here, the black rectangles show the turbines' location.The black dashed lines demonstrate the cyclic boundaries.

Figure 4
depicts contour plots of turbulence intensity at the hub height predicted through LES, the standard k − ε model, and the extended k − ε model.The deficiency of the standard k − ε model in the prediction of the turbulence level in the wake region is visible here.Unlike LES and the extended k − ε model (Figures 4(a) and (c)), the standard k − ε model (Figure 4(b)

Figure 3 :
Figure 3: Distribution of the normalized velocity deficit at hub height in case 1, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε..

Figure 3 :
Figure 3: Distribution of the turbulence intensity at hub height in case 1 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 4 :
Figure 4: The normalized power as a function of turbine row in case 1.Here, normalized power is calculated by dividing the output power of each row by that of the first row through P row /P 1 × 100.The vertical bar indicates the standard deviation of the normalized power obtained from LES.

Figure 5 :
Figure 5: Distribution of the normalized velocity deficit at hub height in case 2 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 4 :
Figure 4: Distribution of the turbulence intensity at hub height in case 1, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 3 :
Figure 3: Distribution of the turbulence intensity at hub height in case 1 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 4 :
Figure 4: The normalized power as a function of turbine row in case 1.Here, normalized power is calculated by dividing the output power of each row by that of the first row through P row /P 1 × 100.The vertical bar indicates the standard deviation of the normalized power obtained from LES.

Figure 5 : 6 Figure 2 :
Figure 5: The normalized power as a function of turbine row in case 1.Here, normalized power is calculated by dividing the output power of each row by that of the first row through P row /P 1 × 100.The vertical bar indicates the standard deviation of the normalized power obtained from LES.

Figure 3 :
Figure 2: Distribution of the normalized velocity deficit at the hub height in case 1, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.The white rectangles show the turbines' location.Figure 3: Distribution of the normalized velocity deficit at hub height in case 1, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε..

Figure 3 :
Figure 3: Distribution of the turbulence intensity at hub height in case 1 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 4 :
Figure 4: The normalized power as a function of turbine row in case 1.Here, normalized power is calculated by dividing the output power of each row by that of the first row through P row /P 1 × 100.The vertical bar indicates the standard deviation of the normalized power obtained from LES.

Figure 5 :
Figure 5: Distribution of the normalized velocity deficit at hub height in case 2 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 4 :
Figure 4: Distribution of the turbulence intensity at hub height in case 1, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 3 :
Figure 3: Distribution of the turbulence intensity at hub height in case 1 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 4 :
Figure 4: The normalized power as a function of turbine row in case 1.Here, normalized power is calculated by dividing the output power of each row by that of the first row through P row /P 1 × 100.The vertical bar indicates the standard deviation of the normalized power obtained from LES.

Figure 5 : 6 Figure 2 :
Figure 5: The normalized power as a function of turbine row in case 1.Here, normalized power is calculated by dividing the output power of each row by that of the first row through P row /P 1 × 100.The vertical bar indicates the standard deviation of the normalized power obtained from LES.

Figure 3 :Figures 9
Figure 3: Distribution of the ratio of eddy viscosity to the inflow value (ν T /ν T 0 ) at the hub height in case 1, predicted by (a 1 ) the standard k − ε and (a 2 ) the extended k − ε.(b) Rotor-averaged ν T /ν T 0 in case 1, predicted by the standard and extended k − ε models.

188 6 Figure 3 :
Figure 3: Distribution of the ratio of eddy viscosity to the inflow value (ν T /ν T 0 ) at the hub height in case 1, predicted by (a 1 ) the standard k − ε and (a 2 ) the extended k − ε.(b) The rotor-averaged ratio of eddy viscosity to the inflow value (ν T,d /ν T 0 ) in case 1, predicted by the standard and extended k − ε models.

Figure 3 :
Figure 3: Distribution of the turbulence intensity at hub height in case 1 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 4 :
Figure 4: The normalized power as a function of turbine row in case 1.Here, normalized power is calculated by dividing the output power of each row by that of the first row through P row /P 1 × 100.The vertical bar indicates the standard deviation of the normalized power obtained from LES.

Figure 5 : 6 Figure 4 :
Figure 5: Distribution of the normalized velocity deficit at hub height in case 2 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 3 :
Figure 3: Distribution of the turbulence intensity at hub height in case 1 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 4 :
Figure 4: The normalized power as a function of turbine row in case 1.Here, normalized power is calculated by dividing the output power of each row by that of the first row through P row /P 1 × 100.The vertical bar indicates the standard deviation of the normalized power obtained from LES.

Figure 5 :
Figure 5: Distribution of the normalized velocity deficit at hub height in case 2 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 5 :
Figure 5: The normalized power as a function of turbine row in case 1.Here, normalized power is calculated by dividing the output power of each row by that of the first row through P row /P 1 × 100.The vertical bar indicates the standard deviation of the normalized power obtained from LES. 161

Figure 6 :
Figure 6: Distribution of the normalized velocity deficit at hub height in case 2, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 6 :
Figure 6: Distribution of the turbulence intensity at hub height in case 2 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.Figure 7: Distribution of the turbulence intensity at hub height in case 2, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 7 :
Figure 6: Distribution of the turbulence intensity at hub height in case 2 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.Figure 7: Distribution of the turbulence intensity at hub height in case 2, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 6 :
Figure 6: Distribution of the normalized velocity deficit at the hub height in case 2, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 6 :
Figure 6: Distribution of the turbulence intensity at hub height in case 2 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 7 :Figure 8 ,Figure 8 :
Figure 7: The normalized power as a function of turbine row in case 2.

Figure 7 :
Figure 7: Distribution of the turbulence intensity at the hub height in case 2, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 6 :
Figure 6: Distribution of the turbulence intensity at hub height in case 2 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 7 :
Figure 7: The normalized power as a function of turbine row in case 2.

Figure 8 :
Figure 8: Distribution of the normalized velocity deficit at hub height in case 3 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 8 :
Figure 8: The normalized power as a function of turbine row in case 2.

Figure 6 :
Figure 6: Distribution of the turbulence intensity at hub height in case 2 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 7 :
Figure 7: The normalized power as a function of turbine row in case 2.

Figure 8 :
Figure 8: Distribution of the normalized velocity deficit at hub height in case 3 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 8 :
Figure 8: The normalized power as a function of turbine row in case 2.

Figure 9 :
Figure 9: Distribution of the normalized velocity deficit at hub height in case 3, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 9 :
Figure 9: Distribution of the turbulence intensity at hub height in case 3 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.Figure 10: Distribution of the turbulence intensity at hub height in case 3, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 10 :
Figure 9: Distribution of the turbulence intensity at hub height in case 3 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.Figure 10: Distribution of the turbulence intensity at hub height in case 3, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 9 :
Figure 9: Distribution of the normalized velocity deficit at the hub height in case 3, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 9 :
Figure 9: Distribution of the turbulence intensity at hub height in case 3 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 10
Figure 10 illustrates the normalized power in case 3 calculated through LES, the standard k − ε model, and our 176

Figure 10 :
Figure 10: The normalized power as a function of turbine row in case 3. 183

Table 1 :
Average normalized power of the waked turbines in cases 1-3 predicted by LES, standard k − ε, and extended k − ε.Relative errors of RANS methods compared to LES are calculated through | LES -RANS |/LES × 100.Error reduction is calculated by subtracting the error of the extended k − ε from that of the standard k − ε.

Figure 10 :
Figure 10: Distribution of the turbulence intensity at the hub height in case 3, predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 9 :
Figure 9: Distribution of the turbulence intensity at hub height in case 3 predicted by (a) LES, (b) standard k − ε, and (c) extended k − ε.

Figure 10
Figure10illustrates the normalized power in case 3 calculated through LES, the standard k − ε model, and our extended version.Owing to the fact that turbine 2 is only affected by a partial wake from turbine 1, its normalized power is significantly higher than the second turbines in cases 1 and 2, and a perfect match between the three models is observed.Focusing on the other waked turbines, the extended k − ε model performs satisfactorily, with a good agreement with the results from LES.

Figure 10 :
Figure 10: The normalized power as a function of turbine row in case 3.

Figure 11 :
Figure 11: The normalized power as a function of turbine row in case 3.

Table 1 :
Average normalized power of the waked turbines in cases 1-3 predicted by LES, standard k − ε, and extended k − ε.Relative errors of RANS methods compared to LES are calculated through | LES -RANS |/LES × 100.Error reduction is calculated by subtracting the error of the extended k − ε from that of the standard k − ε.

2 .
Obtained from the measurements by Chamorro and Porté-Agel [39], to simulate the wind farm, we use the measured values of C ′ T for the different turbine rows as depicted in Figure 11(b).

Figure 11 :
Figure 11: (a) Schematic of a 10 × 3 wind-turbine array.Here, the black rectangles show the turbines' location.The inlet and outlet are located at x/d 0 = −5 and 90, respectively.The black dashed lines show the cyclic boundaries.(b) Disk-based thrust coefficient (C ′ T ) for different turbine rows used for RANS simulation of the 10 × 3 array, obtained from Ref. [41].

Figure 12
Figure12shows the vertical distributions of the normalized streamwise velocity obtained from the standard and extended k − ε models compared against measurement data[39].The profiles are extracted at a distance of 3d 0 behind the turbine rows 1-7 and 10.Focusing on the predictions of the standard k − ε model, one can see that it generally overestimates the normalized streamwise velocity in the area from bottom-to top-tip heights, and this agrees with the earlier observations.In contrast to the standard k − ε model, the agreement between the results of the extended k − ε and measurement data improves consistently from the bottom-tip to the top-tip heights.This observation is also in agreement with previous findings showing that the extended k − ε model predicts more accurate disk-averaged velocities and normalized power of turbines.

Figure 12 :
Figure 12: (a) Schematic of the 10 × 3 wind-turbine array.Here, the black rectangles show the turbines' location.The inlet and outlet are located at x/d 0 = −5 and 90, respectively.The black dashed lines show the cyclic boundaries.(b) Disk-based thrust coefficient (C ′ T ) for different turbine rows used for RANS simulation of the 10 × 3 array, obtained from Ref. [43].

Figure 12 :
Figure 12: Vertical profiles of the normalized velocity at a distance of 3d 0 downstream of rows 1-7 and 10 in the 10 × 3 array of wind turbines, predicted by the standard k − ε and extended k − ε, and compared against data from wind-tunnel measurements.The wind-tunnel measurement profiles (black circles), which relate to the experiments performed by Chamorro and Porté-Agel [39], are extracted from Ref. [41].Horizontal black dashed lines indicate the bottom-tip, hub, and top-tip heights.

Figure 13 :
Figure 13: Vertical profiles of the normalized velocity at a distance of 3d 0 downstream of rows 1-7 and 10 in the 10 × 3 array of wind turbines, predicted by the standard k − ε and extended k − ε, and compared against data from wind-tunnel measurements.The wind-tunnel measurement profiles (black circles), which relate to the experiments performed by Chamorro and Porté-Agel [41], are extracted from Ref. [43].Horizontal black dashed lines indicate the bottom-tip, hub, and top-tip heights.

Table 1 :
Average normalized power of the waked turbines in cases 1-3 predicted by LES, standard k − ε, and extended k − ε.Relative errors of RANS methods compared to LES are calculated through | LES -RANS |/LES × 100.Error reduction is calculated by subtracting the error of the extended k − ε from that of the standard k − ε.