Model-form uncertainty quantification in RANS simulations of wakes and power losses in wind farms

Reynolds-averaged Navier-Stokes (RANS) is one of the most cost-efficient approaches to simulate wind-farm-atmosphere interactions. However, the applicability of RANS-based methods is always limited by the accuracy of turbulence closure models, which introduce various uncertainties into the models. In this study, we estimate model-form uncertainties in RANS simulations of wind farms. For this purpose, we compare different RANS models to a large-eddy simulation (LES). We find that the realizable k-epsilon model is a representative RANS model for predicting the mean velocity, the turbulence intensity, and the power losses within the wind farm. We then investigate the model-form uncertainty associated with this turbulence model by perturbing the Reynolds stress tensor. The focus is placed on perturbing the shape of the tensor represented by its eigenvalues. The results show that the perturbed RANS model successfully estimates the region bounding the LES results for quantities of interest (QoIs). We also discuss the effect of perturbation magnitude on various QoIs.


Introduction
Wind power has been one of the fastest-growing sustainable energy sources in the past few decades and has contributed to the ever-increasing efforts to mitigate greenhouse gas emissions and global warming. Increasing this contribution motivates more attention to larger and more efficient wind farms. This directly leads to computational fluid dynamics (CFD). CFD is a cost-effective approach for modeling and analyzing wind farms, and it can provide detailed information about the wind-turbines-atmosphere interactions. However, it is known that inaccurate prediction of turbulent flow patterns -particularly the wakes formed behind wind turbines -significantly impact CFD's estimates of wind farms' power generation (see, e.g., the review of Refs. [1][2][3]). This justifies efforts to increase the accuracy of the CFD models, on the one hand, and quantifying their uncertainties, on the other hand. In fact, standard CFD approaches are inherently uncertain, and therefore interpretation of any CFD simulation results must be associated with some level of uncertainty quantification (UQ) [4,5].
Different computational approaches have been considered by researchers and engineers, among which eddyresolving simulations -direct numerical simulation (DNS) and large-eddy simulation (LES) -provide a high level of physical fidelity [6]. However, these methods are computationally expensive at high Reynolds numbers, and hence, are more suitable for small-scale problems in fundamental research (see the review of Refs. [7][8][9][10], and references cited therein). Consequently, Reynolds-averaged Navier-Stokes (RANS), a non-scale-resolving tool, is considered the standard method in industrial applications [11,12]. The reduced computational cost of RANS models is associated with the use of RANS closures, which fall short in a variety of physical scenarios, leading to uncertainties and errors in the predicted quantities of interest (QoIs) [4,5]. Specifically, in the wind-energy community, it has been reported in several research publications that standard RANS models fail to accurately predict the wake flow and power loses in wind farms [13][14][15][16], whereas eddy-resolving numerical tools such as LES compare well with laboratory and field-scale experiments [17][18][19][20][21][22][23][24][25][26].
In RANS models, the proportionality between the eddy-viscosity and the mean rate of strain, i.e., Boussinesq's hypothesis, is one of the main sources of uncertainty [4,5]. Since such uncertainty arises due to the modeling assumptions, it is called the model-from (or structural) uncertainly. No matter how much the model constants of these models are tuned according to the problem, this structural uncertainty is indispensable, making it essential to quantify and understand this type of uncertainty in a RANS calculation. One way to quantify the model-form uncertainty in RANS is to decompose the Reynolds stress anisotropy tensor to its amplitude, shape, and orientation represented by turbulent kinetic energy (TKE), eigenvalues, and eigenvectors, respectively. By perturbing any of these values, one can assess the model-form-induced uncertainty in the solution. This method, which was first introduced by Emory et al. [27], has been applied to different flow problems, including wind engineering flows, urban canopy flows, streamlined surface flows [28][29][30][31]. In wind-energy applications, this approach was recently utilized in the RANS simulation of a stand-alone wind-turbine wake [32]. This paper extends the previous work, and we study the model-form uncertainty in wind-farm simulations. Like the previous work that applies the UQ methodology in Ref. [27] to different flow problems, the novelty of this paper does not lie in the methodology itself but rather in the uncertainty estimate. As will be clear in the later sections, many of the results we obtain for a wind farm cannot be directly inferred from the results of a single turbine (as different physics is at play in these two flow configurations).
The specific objective of the present work is to quantify the model-form uncertainty in RANS simulations of wind farms and the associated power losses. The QoIs are the mean velocity deficit, the turbulence intensity, and the turbine power output in wind farms. We will show that the propagation of uncertainties in a wind farm is very similar to that of turbine wakes: they propagate downstream, superimpose, and asymptote shortly after a few rows of turbines. The rest of the paper is organized as follows: section 2 introduces the RANS framework and a description of the numerical setup. A brief description of the LES code is also provided. In section 3, we describe the UQ framework, and in section 4, we show the results. We compare different RANS models to LES and discuss the UQ results. Finally, a summary and concluding remarks are given section 5.

Reynolds-averaged Navier-Stokes (RANS) framework
The RANS framework used in the present study solves the Reynolds-averaged version of continuity and momentum conservation equations for an incompressible flow, where the overbar and prime denote the Reynolds-averaged component and fluctuating component of a variable, respectively, u i is the velocity, ρ is the density, p is the pressure, t is the time, ν is the molecular viscosity,S i j is the mean rate of strain, u i u j is the Reynolds stress, and f i is a body force term representing the wind-turbine effects on the flow field. According to Boussinesq's hypothesis, the deviatoric part of the Reynolds stress tensor is linearly related to the mean rate-of-strain tensor as where k = u i u i /2, δ i j and ν T represent TKE, Kronecker delta, and eddy viscosity, respectively. The computational domain used in RANS simulations is shown in Figure 1, and its size is set to 4400m × 400m × 355m in x, y, and z directions, respectively. In the baseline case, the first turbine is located at a distance of five rotor diameters (D) from the inlet, and the other turbines are located at a 7D distance from the upstream turbine in an aligned configuration. This turbine spacing is quite typical in wind-farm implementations [33], and the aligned configuration promotes wake interactions, whose effect on model uncertainty is a topic of this paper. In addition to the baseline configuration, two additional turbine layouts with staggered turbine arrangement and closer turbine positioning, i.e., 5D will be investigated. The turbines have a diameter (D) of 80m, a hub height (z h ) of 70m, and an axial induction factor of (a) of 0.25. The turbine-induced forces are modeled using the standard non-rotational actuator-disk model described by Calaf et al. [34] using the wind velocity at the rotor plane and the disk-based thrust coefficient, , which is set to 4/3.
The inlet and top boundary conditions are defined by prescribing velocity, turbulence kinetic energy, and turbulence dissipation rate (or specific dissipation rate), corresponding to a logarithmic boundary-layer flow, together with a zero-gradient condition for the pressure [35]. The outflow boundary condition is a pressure outlet with zero gradients for the velocity and other turbulence quantities. A cyclic boundary condition is applied in the lateral direction for all the variables. A rough wall boundary condition with wall functions for the turbulence quantities is imposed at the bottom wall [16]. The inlet's velocity at the hub height (ū h ) is 8m/s, and the turbulence intensity (I = √ 2k/3/ū h ) at the same height is 5.8%. The aerodynamics surface roughness (z 0 ) is determined for each turbulence model to satisfy the desired conditions at the hub height. The flow's governing equations are discretized and solved using the OpenFOAM open-source software, based on the finite volume approach and the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm [36]. Different turbulence models considered in this study are described in the next section. A grid convergence study is conducted, and the results are shown in Appendix A.

Large-eddy simulation (LES) framework
Here, an in-house pseudo-spectral code is employed for the LESs, which has been well validated and used in earlier wind-energy research publications (e.g., Refs. [37][38][39][40]). The size of the computational domain is 4800m × 800m × 355m in x, y, and z directions, respectively, and it is discretized uniformly into 480 × 160 × 72 grid points. The LES domain in the x-direction is slightly larger than the RANS setup due to the presence of a fringe zone, which is used to adjust the flow from the downstream wake state to an undisturbed inflow condition [41]. There are two columns of six turbines with the same streamwise spacing as in the RANS simulations and with the lateral spacing of 5D. The width of the RANS domain is half of the LES domain, and LES results are averaged over the two columns to be compared to RANS. For the sake of brevity, further details about the LES code have been excluded, as well as its equations here. A more detailed description of the LES code can be found in Refs. [42,43].

Uncertainty quantification (UQ)
The primary purpose of this study is to investigate the model-form uncertainty in our RANS simulations. In the following, we describe the UQ methodology. We represent the Reynolds stress tensor in terms of its amplitude, shape, and orientation [27,44]. The Reynolds stress anisotropy tensor (b i j ), which is the deviatoric part of the Reynolds stress normalized by the TKE, can be written as [45] The eigendecomposition of the normalized anisotropy tensor yields (see, e.g., Ref. [46]) Here, Λ mn is a diagonal tensor containing the eigenvalues of the anisotropy tensor in a decsending order (i.e., λ 1 ≥ λ 2 ≥ λ 3 ), and v i j is a tensor containing the orthogonal eigenvectors of the anisotropy tensor. Note that due to the normalization constraints, the trace of the diagonal eigenvalues tensor is equal to zero, i.e., λ 1 + λ 2 + λ 3 = 0 [45]. Using Equation 3 and Equation 4, the Reynolds stress tensor can be written as Applying perturbations to the amplitude (k), shape (Λ), and orientation (v) of this tensor, we have the following perturbed Reynolds stress tensor where k * , v * im and Λ * mn are the perturbed TKE, eigenvectors and eigenvalues, respectively. Adding the difference between the perturbed and the original Reynolds stress tensors (i.e., ∆R i j = u i u j * − u i u j ) to the RANS equations gives rise to the perturbed RANS equations ∆R i j is the model-form uncertainty. Following Refs. [30,32], we implement Equation 7 in the open-source solver OpenFOAM. Here, we focus on perturbing the eigenvalues of the anisotropy tensor, which can be represented on a Barycentric map. Having a triangle, each of its corners represents a limiting state of one-componentâ 1 , twocomponentâ 2 , and three-component (or isotropic)â 3 turbulence. One can show that any Reynolds stress can be represented on the Barycentric triangle map by a linear combination as are the coordinates of the Barycentric triangle's corners [27]. The shape of the Reynolds stress tensor can be perturbed towards one-, two-, and three-component limiting states as where δ is the amount of perturbation from 0 to 1, and Λ c could be one of the limiting states defined as It worth mentioning that by assigning an RGB color code to each of the points inside the Barycentric triangle, the Reynolds stress tensor can be visualized at any arbitrary location. The color interpretation used in this study is shown in Figure 2. The following equation calculates the RGB code for any point in the Barycentric triangle map as where C 1 = (λ 1 − λ 2 ) + 0.65, C 2 = 2(λ 2 − λ 3 ) + 0.65, and C 3 = (3λ 3 + 1) + 0.65 [47]. Here, the coefficients of 0.65 and 5 are for optimal visualization and better combination of colors.

Evaluation of different RANS models
Among different closure models that can be used to determine the Reynolds stress term, the two-equation eddyviscosity models are the most commonly used by the wind-energy community, and this study focuses on two-equation models. In this study, we assess the performance of five different closure models: standard k − ε [49,50], fp k − ε [51], RNG k − ε [52], realizable k − ε [53], and SST k − ω [54]. Note that a detailed comparison of RANS models is not the purpose of the present work as it has been addressed in the previous studies (e.g., Refs. [55][56][57][58][59], among others). The purpose here is to select a baseline RANS model for the UQ study. The details of the turbulence models (their equations and model constants) are omitted for the sake of brevity, and the reader may refer to the above-mentioned references for further information. In the following, the results of different RANS simulations are presented and compared with the LES data. The focus is placed on the prediction of normalized velocity deficit, defined as ∆ū/ū h = (ū in −ū)/ū h , turbulence intensity, and power outputs in the wind farm. Here subscripts in denotes inflow.
The simulated normalized velocity deficit is shown in Figure 3 in a two-dimensional horizontal plane at the turbine hub height. The normalized velocity deficit, averaged over the rotor area, is also plotted as a function of the downwind distance in Figure 4. Both figures indicate that the RANS models considered here underestimate the wake velocity deficit behind the most upstream wind turbine. The SST k-ω model has a comparably better performance behind the first turbine, but its prediction deviates from the LES for the downstream waked wind turbines. The largest deviation, especially behind the first turbine, is observed for the standard k-ε, which can be attributed to the overestimation of the eddy viscosity immediately behind the turbine [13,51]. We see from the two figures that, for the particular case considered here, the RNG k-ε, the fp k-ε, and the realizable k-ε models show relatively good agreement with the LES data in the near wake and far wake for all downstream turbines. Figure 5 shows the contours of the turbulence intensity at turbine hub height. The turbulence intensity, averaged over the rotor area, is also plotted as a function of the downwind distance in Figure 6. We see that the standard k-ε does not do well in term of estimating the turbulence intensity in the wake of the first turbine. In particular, unlike the other RANS models, it fails to capture the double peak structure of the turbulence intensity behind the first turbine. The fp k-ε performs well behind the first turbine, but it significantly overpredicts the turbulence intensity downstream of the other turbines. It can also be seen that the SST k-ω and RNG k-ε models underestimate turbulence intensity for downstream waked turbines. Overall, we see that compared with other models, the Realizable k-ε shows a better agreement with the LES data in estimating the turbulence intensity.
A comparison of different RANS models in predicting wind-farm power output is provided in Figure 7. The results are normalized with the power of the most upstream wind turbine. Consistent with the results shown in Figure 3 and Figure 4, it can be seen that the standard k-ε and the fp k-ε overestimates the power output for the waked turbines. SST k-ω model provides a good estimation for the power output of the first downstream turbine but underestimates the power for the rest of the downstream turbines. RNG k-ε model shows a similar trend as SST k-ω model. The realizable k-ε shows a better agreement with the LES in predicting the power output, although its prediction of the efficiency of the first downstream wind turbine is not entirely accurate.     Overall, we find the realizable k-ε model performs well and is a representative RANS model. Hence, it is selected as the baseline model for the UQ study.

Uncertainty quantification (UQ) analyses
Before we proceed with detailed UQ analysis, we first examine the Reynolds stress anisotropy in both RANS and LES. First, Figure 8 shows the RGB colormap of the Reynolds stress anisotropy tensor at the horizontal plane through the hub height for LES and the baseline RANS model following Equation 12. As can be seen, near the two side-tip positions, the Reynolds stress obtained from LES tends to be more anisotropic (towards the one-component turbulence) compared to the RANS results. This is expected as the RANS eddy-viscosity models underestimate the level of anisotropy in the regions with a high level of shear. It is also observed that, in the core of the wake, the baseline RANS model underestimates the level of isotropy compared to the LES data.
To have a more quantitative insight on the results presented in Figure 8, the structure or the Reynolds stress anisotropy on the Barycentric map is explored in Figure 9 at 5D downstream of each individual turbine and at the hub height level, for both LES and RANS models. The results show a fairly similar pattern for all six locations. In all  The results shown in Figure 9 indicate that a moderate perturbation of the Reynolds stress tensor towards onecomponent and three-component turbulence will make the RANS predicted Reynolds stress more similar to the LES under this particular metric, and that a perturbation towards the two-component extreme will make the RANSpredicted Reynolds stress less like the LES. However, a close agreement in the Barycentric map does not guarantee whatsoever a close agreement for other QoIs. That is, one cannot infer from Figure 9 how a perturbation towards the one-, two-, three-component extremes will affect RANS's prediction of QoIs like the velocity deficit, the turbulence intensity, and the power loss. As will be shown later, perturbation of the RANS model toward the two-component turbulence leads to a prediction that lies within the ones obtained from the other two limiting states of turbulence. Here, three different values for the magnitude of perturbation δ are considered, i.e, 0.25, 0.5, and 1. The δ = 0.5 is adopted based on a priori studies [32], which is also in good agreement with the results presented here. In addition to this value, one extremely conservative (δ = 1) and one less conservative value (δ = 0.25) are also considered to analyze the effect of injected perturbation amount on the QoIs. Here, we follow the previous work and use a uniform δ. We can get a non-uniform δ by make it a random number, which may not be very meaningful. A more meaningful way of getting a non-uniform δ is to make it a function of the local flow condition [60]. This, however, is usually only attempted when calibrating a model, in which case, the function is obtained through fitting. As model calibration falls out of the scope of this work, we follow the previous study and use a uniform δ. Figure 10 represents normalized velocity deficit contours for δ = 0.5 at the horizontal plane through the hub height. As can be seen, perturbation towards one-component turbulence increases the velocity deficit, leading to a more extended wake. In contrast, perturbation toward the three-component (isotropic) state enhances the mixing and leads to a faster wake recovery. It is also observed that perturbation towards the two-component turbulence yields a prediction that lies within the ones obtained from the one-component and three-component states of turbulence. Figure 11 shows the effect of Reynolds stress perturbation towards three limiting states on the turbulence intensity for δ = 0.5. As expected, perturbation towardsâ 1 increases the shear and, consequently, the turbulence intensity behind the turbines. In contrast, the perturbation toward the three-component state leads to a reduction in the turbulence level. These results are consistent with the ones obtained for the wake behind a stand-alone wind turbine [32]. Similar to the results for normalized velocity deficit, model perturbation towards the two-component turbulence lies between the predictions obtained from one-component and three-component states.
The evolution of the normalized velocity deficit, averaged over the rotor area, is shown in Figure 12 for different amounts of perturbations. As can be observed, perturbation toward one-component turbulence increases the velocity deficit; perturbation toward the three-component state decreases the velocity deficit; and perturbation toward the twocomponent turbulence lies within the other two states. It is also evident that increasing the δ value from 0.25 to 0.5 and subsequently to 1.0 increases the uncertainty bands in the RANS prediction. It can be seen that, except in the near    wake region behind the first turbine, the LES results are well covered in the band between the one-component and three-component curves in Figure 12. Figure 13 shows the variation of turbulence intensity, averaged over the rotor, with downstream distance for different values of δ. Consistent with the previous results, perturbation towards one-component and three-component turbulence provides the upper and lower bounds in predicting the turbulence intensity behind the turbines, respectively. Again, it is observed that the higher the amount of perturbation, the wider the bounds. It is also found that all three different perturbation bounds capture the LES data fairly well, except in the vicinity of the turbines, in the near wake region. This might be related to the fact that the baseline model shows poor performance in capturing the magnitude of TKE immediately behind the turbines.
For a more detailed analysis of the effect of perturbation on the QoIs, the lateral profiles of the velocity deficit and turbulence intensity at 5D downwind of each individual turbine are shown in Figure 14 for δ = 0.5. As can be seen, the velocity deficit and turbulence intensity obtained from LES are covered relatively well by the bounds obtained from one-and three-component perturbations. Consistent with the previous results, perturbation towards one-component and three-component states respectively increases and decreases the velocity deficit in the rotor region. Also, as can be seen here, perturbation towards the one-component state reduces the turbulence mixing, yielding stronger the double-peak structure for the turbulence intensity, and a less smeared wake deficit.
The impact of model-form uncertainty on the predicted power output of the wind farm is shown in Figure 15.  Consistent with the previous results, perturbation towards the one-component turbulence leads to an increase in the velocity deficit and, consequently, a decrease in the power output of the downstream turbines. On the other hand, perturbation toward three-component state enhances the wake recovery and, thus, increases the power output of the waked turbines. It can be seen that when the amount of perturbation value increases, the uncertainty bounds widen. It is also found that the model perturbation using the suggested value of δ = 0.5 leads to accurate predictions of the power output.
To shed light on the impact of wake interactions on uncertainty bounds, the normalized velocity deficit and turbulence intensity bounds for baseline and single-turbine cases are compared in Figure 16. As can be seen in this figure, for a single turbine, the uncertainty bound is larger in the region associated with the higher turbulence intensity and decreases with the downstream distance as the wake recovers in the far wake region. However, in the case of the wind farm, the region behind each individual turbine is associated with a high level of shear and turbulence and, consequently, a larger uncertainty bound.
Finally, in order to further explore the model-form uncertainty on the power losses of wind farms, two additional layouts for the wind farm are considered. In the first configuration, turbines are aligned similar to the baseline case, but the streamwise spacing between them is 5D representing a case with stronger wake interactions. In the second Fig. 16: Normalized velocity deficit and turbulence intensity, averaged over rotor area, for LES, Baseline RANS, and perturbed RANS towards one-, two-and three-component turbulence with δ = 0.5. From top to the bottom: normalized velocity deficit through the wind farm, normalized velocity deficit through the stand-alone turbine, turbulence intensity through the wind farm, and turbulence intensity through the stand-alone turbine.
configuration, turbine distances in the streamwise direction are 7D, but the turbine rows 2, 4, and 6 are shifted by 1D to the right side with respect to the most upstream one. The reason for selecting the latter case is to show the performance of the proposed methodology to the case with partial wake interactions. Since the results for velocity deficit and turbulence intensity show a similar trend as the baseline case, and for the sake of brevity, we only show the results related to the power losses. As can be seen in Figure 17, similar to the baseline configuration, perturbation towards one-and three-component turbulence, with δ = 0.5, can bound the LES results for both cases relatively well. These results show the robustness of the method for quantifying the model-form uncertainty in the prediction of power losses of wind farms under both full-and partial-wake conditions.

Summary and concluding remarks
This study aims to quantify the model-form uncertainty in the RANS simulation of wakes and power losses in wind farms. The performance of different RANS models is first tested with the LES data in the prediction of mean velocity deficit, turbulence intensity, and power output from wind turbines. Different RANS models considered in the study are as follow: the standard k-ε [49,50], fp k-ε [61], RNG k-ε [52], realizable k-ε [53], and SST k-ω [54]. It is found that, the realizable k-ε model shows a relatively good agreement with the LES data for the above-mentioned quantities, and it is selected as the baseline model for the UQ study.
In order to quantify the model-form uncertainty in the baseline RANS closure, the state of turbulence is perturbed towards one of the three limiting states; one-component, two-component, and three-component (isotropic) turbulence. In this study, three different values for perturbation amount δ= 0.25 (the least conservative coefficient), 0.5 (obtained from a priori studies [32]), and 1 (extremely conservative approach) are considered. Simulation results show that perturbation towards the three-component turbulence leads to a faster wake recovery and reduces the turbulence intensity within the wind farm. Consequently, the power output of the downstream turbines increases when the structure of turbulence is perturbed toward the three-component state. On the other hand, perturbing the turbulence state towards the one-component turbulence leads to an opposite trend: increasing the wake velocity deficit and turbulence intensity and reducing the power output of the waked turbines. Perturbation towards the two-component state does not significantly change the model predictions compared to the baseline model, as it lies within the results corresponding to the one-component and two-component perturbations. It is also found that the value of δ = 0.5 is a suitable choice for creating an acceptable bound covering the LES results. This is further confirmed through two additional wind-farm configurations with staggered turbine arrangement and closer turbine positioning.
In the present study, the focus is placed on perturbation of the shape of the Reynolds stress tensor represented by its eigenvalues. Future work can consider the application of Global UQ [62] and data-driven machine learning [63] to improve the accuracy of the RANS-predicted Reynolds stress tensor in wake-flow simulations. In addition, the use of non-uniform perturbation magnitude as well as non-uniform regions for perturbations can be further studied (see, e.g., Refs. [64,65]). Future research is also required to address the effect of model-form uncertainties in RANS simulations of wakes and power losses under thermally-stratified conditions and/or over complex terrain.